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M. A. BIOT 
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Abstract. Closed form expressions are developed for the output of a frequency 
modulation receiver for an arbitrary number of superposed input signals. This corresponds 
to problems of interference or disturbance due to scatter and multiple reflexions. It is 
also shown how the Fourier components of the output may be evaluated by methods 
more direct than the usual Fourier analysis. 

Introduction. A frequency modulation receiver is essentially a non-linear device. 
The input signal is usually fed into two non-linear filters, first into an amplitude limiter 
which reduces the signal to a constant amplitude, then into a discriminator whose 
output is a rectified signal with an amplitude proportional to the frequency deviation 
from the carrier frequency. Sometimes this output of the discriminator is processed 
through a linear filter which eliminates all but a few frequency components of the 
modulation. 

Because of the non-linearity of the system, special methods must be devised to 
evaluate the output due to the superposition of input signals. The procedure presented 
here yields closed form expressions for the output when an arbitrary number of signals 
or a continuous distribution of them are superposed. The results may be used to predict 
interference effects or the disturbance due to multiple reflexions or scattering of the 
main signal. The method makes use of the concept of instantaneous frequency. The 
limitations of this concept in analyzing the behavior of frequency modulation circuits 
was discussed extensively by Carson and Fry‘ and Van der Pol’. 

The general theory is developed in Section 1 for the case of an arbitrary input repre- 
sented by a continuous spectrum. This is applied in Section 2 to the case of an arbitrary 
number of signals with a single modulation frequency. It is also indicated how the 
method applies when there is more than one modulation frequency. Section 3 deals 
with the Fourier analysis of the output. Because of the fact that the expression for 
the output is in the form of the quotient of two Fourier series, methods more direct 
than the usual Fourier analysis are applicable. It is shown that the Fourier coefficients 


*Received March 8, 1956. This report covers work performed under Contract DA-30-115-ORD-47, 
Ordnance Project No. TUI-2020, Dept. of Army Project No. 516-05-002 with the Rochester Ordnance 
District, Dept. of the Army. 

1J. R. Carson and T. C. Fry, Variable frequency electric circuit theory with application to the theory 
of frequency modulation, Bell Syst. Tech. Journal 26, 513-540 (1937). 

2B. Van der Pol, Fundamental principles of frequency modulation, Journal I.E.E. 111, 153-158 (1946). 








2 M. A. BIOT [Vol. XV, No. ] 


may be evaluated directly by use of the theorem of residues of analytic functions. The 
method yields directly the Fourier components of the output in both amplitude and 
phase for any number of superposed input signal. In particular, the input may be repre- 
sented by the superposition of its spectral components. This has an important bearing 
on scatter problems since it is then reduced to the calculation of the scatter for each 
monochromatic component. The procedure has been applied to a number of practical 
cases and found to be quite satisfactory from the standpoint of simplicity and accuracy. 
These applications along with some further simplifications will be presented in a sub- 
sequent paper. 

1. General expressions for the output signal of an FM receiver. We consider an 
input signal which is both amplitude and frequency modulated, 


E(t) = 31()[exp (i¢) + exp (—7¢)] (1.1) 
with 
> = wt + g(t), 


where w, represents the carrier frequency, g(t) the frequency modulation and /(t) the 
amplitude modulation. In the receiver the amplitude is first reduced to a constant 
value in a limiter circuit. It is then processed through a discriminator. This is usually 
made of circuits slightly off resonance with the carrier frequency, such that the output 
is a rectified voltage proportional to the frequency deviation of the signal from the 
carrier frequency. This frequency deviation being d¢/dt — w, , the ouput of the receiver 


is then given by 
dp ) , dy 
- a x=. Iz 
x(% a ale” CH 
Sometimes this output signal is passed through a linear filter so as to extract frequency 
components in a narrow band. 
In practice, the signal is not always given in the form (1.1) so that one cannot use 
(1.2) to evaluate the receiver output. In particular, we are interested in the calculation 
of the receiver output when the incoming signal is given by its spectrum 


E(t) = im Glw)e'** dv. (1.3) 


We establish a mathematical processing of the expression (1.3) by which it is possible 
to evaluate a quantity proportional to the rate of change of the phase angle ¢, hence 
also proportional to the receiver output. To do this we apply to the signal a linear operator 
which consists in replacing its spectrum G(w) by 


G@w)(1 + KQ), (1.4) 
where 
2 = |w|] — w,* 


is the frequency deviation of the spectral component from the carrier frequency. The 
effect of this linear operation on the signal can be readily evaluated in the form (1.3) 


*If we wanted to restrict ourselves to analytic functions, the same purpose could be accomplished 
by putting 2 = (w? — «? 


Qu. . 


) 
e// 
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by multiplying G(w) by 1 + KQ in the integral. We shall now introduce the basic as- 
sumption of the method, namely that this operation is approximately equivalent to 
multiplying the signal by the factor 


= 1 +K¥. (1.5) 


This assumption is essentially the same as that upon which is based the design of a 
discriminator circuit, namely that the response of the circuit to the instantaneous 
amplitude and frequency of the signal is the same as in a steady state. The assumption 
will, of course, apply if the frequencies at which the amplitude, 7, and phase angle, ¢, 
vary are very much smaller than the carrier frequency, », . 

If we consider the signal in the form (1.1), it becomes 


E,(t) = }AI(t)[exp (i) + exp (—74)]. (1.6) 
Squaring this quantity we obtain 
EX(t) = 3A°T’ + tA°l’[exp (2i¢) + exp (—279)]. (1.7) 


We notice that the first term represents the low frequency components, while the second 
term represents the high frequency components. In practice, these components are widely 
separated. We may write 


LEV) = 3A°T’, (1.8) 


where the symbol.£ signifies “low frequency part of ---”. Similarly, if we square the 
original signal E(t) we may write 


LE*(t) = 41’. (1.9) 
Hence, 

£Ex(t) = A’? = IK (22). 

LE) ~ A*=1+ ox & E+ K dt (1.10) 


The linear term in K in the expression is proportional to the output signal of the receiver. 
We may write this output signal as 


~cxw¢w__K_|d om 
M() = K% = atal& oxi | . (1.11) 


This expression will now be evaluated in terms of the representation (1.3) of the input 
by means of a spectrum. The square of the signal is the double integral. 





EX) = [- i GOGw)e"** dé de. (1.12) 


The integral is extended to the infinite plane. We introduce the following change of 
variables 


E=3n+9, w= in- $) (1.13) 


E%(t) = tf [- a? . a(2 5 Ee dn dt. (1.14) 


and derive 
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The spectrum of 2E7(¢) is therefore 





R(n) = / o(? as a(2 = ‘) dt. (1.15) 
In evaluating this expression we take into account the fact that G(w) is small except 
in the vicinity of w = + w, . Hence, the contribution to the integral will be only in 


the vicinity of the four points (see Fig. 1). 





“ a 
cy | \ A 
Fat 
Fia. 1 
fF; 9 = 2w, ie = 0 
P,: n=-W. ¢=0 
P;: »=0 o = 2, 
Py: 7 =0 f= —2w, 


The low frequency components of E’(t) are therefore given by integrating (1.14) in the 
vicinity of points P; and P, 


+2 1 — ¢ t Y > 
LE(t) = = I ott fe) dS, (1.16) 
PsP. . - 


— ~ 


the integral being evaluated over elements of area dS in the vicinity of P; and P, . 
Similarly, the low frequency components of E;(t) are given by 











£Ei(t) = ‘.. 1+ K(/w|—o,)J[1 + K(/é|—- w)i@(2t)o( 2H) dS. (1.17) 

Hence, 

1 d 2 

2 Kl a rx, (1.18) 
nfo, (Apt leet laces) (eas 

From (1.16), (1.18), and (1.11) we derive the receiver output 


. -1 
yf 7 $ \nli 7 i Y 
{ff ee (2+ )e (45 ) " as} 
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We have thus expressed this output in terms of the input spectrum G. It is seen that, 
as required by the problem, this expression is independent of the input amplitude. 
Expression (1.19) may be put in a somewhat different form by introducing 


(EH tae 


mf (Lett wolfe) a. 


2 





F(n) 


Il 





J—a@ - 


The equivalence of these two expressions is easily seen if we replace the variable ¢ by —¢. 
This function constitutes the spectrum of the numerator of expression (1.19). With the 
spectrum R(n) of the denominator, as defined by (1.15), we may write 


M(t) = K | F(n)e'™ anf R(n)e"™ any . (1.21) 


The integrands in the expressions for F(y) and R(7) are different from zero only in the 
vicinity of £ = + 2w, and the variable 7 is restricted to the vicinity of the origin. The 
range of integration — « < 7 < + eis that of the low frequency portion of the spectral 
functions F(n) and R(m). According to the footnote remark (*), we could also write 


for F(n) 





F(y) = [ 3 E igo = ot lo(2 + Sa(2 5 ‘) dt. (1.22) 


=. “% 


This expression will be approximately equal to (1.20). 

2. Application to the superposition of sinusoidally modulated signals. In certain 
cases the input signal is made up of the superposition of signals whose frequency modu- 
lation is sinusoidal with a common carrier frequency, but with a different phase for 


each modulation. The signal is then expressed as 


2E(t) = >> R; exp & + 19; + i() sin (w,t + | 
i ™ (2.1) 


+ > R* exp | — it Spe 19; bad (4) sin (wit + |, 
] 1 


where R; and R* are complex conjugates. Such a signal occurs, for instance, in the 
case of multiple reflection or scatter. In this case the phase differences are 


9; = —wl; , (2 2) 
Vi _ —w,l; ’ 
where ¢; is the time lag for arrival of the jth component in the receiver. In order to 
apply the results of the previous section, we must represent the signal by its spectrum. 
We make use of the identity 


exp (78 sin 6) = > J,(8) exp (nié), (2.3) 


n=—@ 


where J, is the Bessel function of the first kind of order n. By putting 8 equal to the 


modulation index 
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s= (2.4) 
w; 
and 
B, = >); R; exp (iy; + iny,) (2.5) 


Bt 


Dui Rt exp (—ig; — iny) 


we may write the signal as 


2E(t) = >> B,J, exp (iw,t + inw,t) ++ D> BtJ, exp (—iw.t — inw,t). (2.6) 


n=—© n=—@ 


This latter expression constitutes the expansion of the signal into a discrete spectrum of 
equidistant frequencies. The integration (1.3) is here replaced by a summation. The 
low frequency components of the square of the signal corresponding to expressions 
(1.16) are given by 

42E°(t) = >) B,BtJ? 


n=—@ 


+ exp (iw,t) >> B,B* JnJn-1 + exp (—tw,/) >, BBt.Jdder (2.7) 


+ exp (2iw,t) >> BBY 2JaJn-2 + exp (—2iw,t) > BBY eJacs 


ga—@ n=—@ 


+ ete. 


We note that 


DL BBiIJarr = LD Bar Btn 


+o +@ 9 
Y BBiaJwer = D BysBtJsoJs — 
etc. 
and put 
C, - > BBs id ad a2 
Tor (2.9) 
ce saa + BoB iW ad as ° 
Then (2.7) may be written 
42E*(t) = Co + Do [Cre***! + Cte **"*]. (2.10) 
k=1 


We must also evaluate Z(t). This is obtained by multiplying by 1 + KQ each frequency 


component in the expansion (2.6) of E(t). In this case 
—— i. 
oe ie) mt (2.11) 


w = +(w, + mw). 
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We note that in practice the terms in the series (2.6) are vanishingly small for | n| > 8 
because, in that case, J, & 0. It is therefore legitimate to write 
Q = nw, . (2.12) 
Putting 
A, = 1+ Kae, (2.13) 


we find 


2E,(t) = >> A,Bul, exp (iwet + inw,t) + D> A,BtJ, exp (—iw.t — inal). (2.14) 


n=—@ n=—@ 


Proceeding as before 


42E%(t) = Fy + > [F, exp (tkw,t) + F% exp (—ikw,)] (2.15) 
k=1 
with 
F, = >> A,A,-2B, Btu id ad 2-2 
=i (2.16) 
Ft - pi A,A,-1B*By-id ad n-2 . 


In order to obtain d/dK &Ej(t) we consider the factor A,A,-, which is the only one to 
contain K. We have 


A,A,-. = [1 — nKw,][1 + (n — k)Ka,] (2.17) 

and 
| A,A | = (2n — k)w (2.18) 

dK 4in+in-k nite ~ ie ° 
Hence, 
sl 2 ox) | = w,Hy + w, > [H, exp (ikw,t) + Ht exp (—ikw,t)] (2.19) 
K-0 k=l 
with 


ps (2n aie k)B, Bt id nd ns 
oi (2.20) 


+0 


Ht = >> (2n — KBtB, WJ aJa-: - 


a=—@ 


A, 


The output signal due to the input E(é) given by (2.1) is 


— =F | 
M(t) = SEED E LE%(t) a 


k 


Ko, 
2 


-—@ 





Co + > [Cy exp (shut) + Ct exp (—iw:kd)] 
k=-—@ 
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The same method may be applied if the superposed signals are not modulated by a 


single sinusoidal component. Consider, for instance, a signal component E;(¢t) con- 


taining two simultaneous modulation frequencies w, and w, . 
E,(t) = exp [iw.t + ig, + ig.] + exp [—iw.t — ig, — tg] (2.22) 


with 


Il 


¥1 


B, sin we (2 23) 


Qo Bo SIN wel. 
Using the identity (2.3) the spectrum of E,(¢) is obtained by writing 
E(t) = exp iw.t exp ig, exp t¢2 


exp 1, (3 > J (81) exp (tnw,t |: : J,(82) exp (tnw. 0 | 


n=+o 


(2.24) 


Il 


Performing the multiplication yields a discrete spectrum. However, this time the 


frequency intervals are not equal. From (2.24) we derive the spectrum due to the super- 
position of signals of the type (2.22) with individual time lags and amplitude factors 


as in (2.1) | 
E(t) = D RE — t,). (2.25) 


Proceeding as above, the spectrum may be used to evaluate the receiver output due to 
this superposition. 

3. Fourier analysis of the receiver output. We consider the case of the superposition 
of input signals modulated by the same modulation frequency w, . We have shown that 


the output is given by expression (2.21) which is the quotient of two Fourier series. It 
which we may expand into a Fourier 


is, itself, a periodic function of frequency «, , 
series. We omit the constant factor in expression (2.21) and write, putting w,f = r 


H, + >i [H, exp (ikr) + H* exp (—tk7)] 
Mi) « —~—221.. 
Co + > [C, exp (tkr) + C% exp (—ik7z)] (3.1) 
k=1 


M, + > [M, exp (iks) + Mt exp (—iks)]. 
k=1 


If the receiver output is filtered through a lowpass filter so that only the fundamental 


component of the Fourier series of M is observed, we must evaluate the coefficients 


M, and M*%. We have 


a2r 


2nM, = | M(r) exp (—?tr) dr 
; (3.2) 


a2 


2rM* = | M(r) exp (tr) dr. 
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The integrals may be evaluated by a method which takes advantage of the particular 
form of the function M/(é) in the present case. Consider the second integral and change 
the variable of integration 7 to a complex variable p 


p = exp (i7) (3.3) 
and express M(r) in terms of p 


H, + > [H.ip* + Htp™“*] 








M(r) = a (3.4) 
Co + > [Cip* + Ctp™] 
=1 
The value of M% is then given by the contour integral on the unit circle 
, p Hot 2 (hip + Ht] 
M* = Oni at dp. (3.5) 


Co+ D> [Cixp* + Ctp™] 
k=1 


The value of this integral is equal to the sum of the residues and depends on the poles 
contained within the unit circle. There are multiple poles at p = O and poles corre- 
sponding to the roots of the denominator. These are the roots of the equation 


Cot > [Cip* + Ctp™*] = 0. (3.6) 
If 
p, = 7, exp (76) (3.7) 


is a root of this equation, then 
Co + } [Curt exp (k0,) + Ctrz* exp (—ik6,)] = 0. (3.8) 
The conjugate of this expression must also be zero, hence, 
Co + : [Ctr; exp (—7k0,) + C.ri* exp (iké,)] = 0. (3.9) 
But this expresses that 
P2 = - exp (76) (3.10) 


is also a root of equation (3.6). Therefore, the roots are grouped in pairs of the same 
argument @, and such that the product of their moduli is unity. Half the roots will be 
inside the unit circle and to each of these roots corresponds an outside root on the same 
radius from the origin. If there are roots on the circle, the above conclusion does not 
hold. However, since the expression (3.6) on the circle represents the square of the 
amplitude of the input signal J*(¢), this can only happen if the input vanishes. As an 
example, let us consider the case when the denominator contains only C, and C, . This 
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will often be the case in practical applications when higher order terms are negligible. 
Denote by N(p) the numerator of the integrand in (3.5). This integral becomes 


1 N(p) dp 
ae g —. 

. 2r1 Co + Cip + Ctp™ (3.11) 
The roots of the denominator are 


—C, + (C2 — 4€,C%)"” 
—r = ——- 





(3.12) 


—C, — (C5 — 4€,C%)'” 
P2 20, Z 





We shall assume the radical is real, in which case the roots are not on the unit circle. 
Moreover, if p, is inside the unit circle, the other root p, is outside. 
The integrand of (3.11) is 


Hop + > (H.p**' + Htp**'] 
k=1 





— — (3.13) 

Ciyp + Cp t+ Ct Ci\(p — pi)(p — pro) wile 
There is a residue due to the root p, and residues due to the terms Htp-**’ for k > 1 
corresponding to poles of order k — 1 at the origin. The sum of residues inside the unit 
circle is 





1N(p,) = H* | = ( 1 )] ; 
Mt = eS | ang (a ee . (3.14 
Ci(pi — pz) bd 2 (k — 2)! Ldp** \C\p’ + Cop + C%4/ J, <0 ) 


This expression gives the phase and amplitude of the fundamental Fourier component 
of the output signal. If the denominator contains more terms than assumed here, we 
must solve a complex algebraic equation of higher degree and similarly evaluate the 
residue for these roots inside the circle. If the numerator N(p) contains only one Fourier 


component 
N(p) = H, + Hip + H%p™' (3.15) 


then (3.14) reduces to the very simple form 


* 7 2 
_ Hi + Hop: + Asp. (3.16) 


* 
M* he ee “cu 


Instead of using the theory of residues, another method of computing the coefficients 
M, of the Fourier expansion (3.1) is to write 


H, + ss [H.p* + Htp*] 





= = M,+ >> [M.p* + Mtp“] (3.17) 
- k= 
Cot+ DY [Cup + Ctp“] ; 

k=1 


considering M, as undetermined, then to multiply both sides of this equation by the 
denominator and equating coefficients of the same power of p. 


ll 


ON THE SOLUTION OF A DIFFERENTIAL EQUATION WITH NONLINEARITY 
APPEARING IN THE SECOND DERIVATIVE OF COMBINED LINEAR AND 
CUBIC TERMS* 


BY 
CHI-NENG SHEN 
Thayer School of Engineering, Dartmouth College 


1. Introduction. One commonly analyzed application of automatic controls to 
processes is shown schematically in Fig. 1. Here a linear feedback controller is used to 
control a two-capacity process. The process involves two linear capacities or tanks in 
this case, and two linear resistances. The analysis of this process is treated in detail in 
the literature and in the texts on automatic control. The solution is simple and straight- 
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Fig. 1. Linear two-capacity process with controls. 
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Fic. 2. Single nonlinear capacity of two-capacity process. 
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forward. If, however, nonlinear elements are introduced, the problem becomes more 
complex and cannot be solved in the conventional manner. 

If tanks of the form shown in Fig. 2 are used, the capacity is nonlinear. The capacity 
of a parabolic tank may be expressed by an equation of the form: 


Volume = C(x + az® + yx"). (1) 
If y is small, this equation reduces to: 
Volume = C(x + az’). (2) 












__ ,, _PLanel oF SYMMETRY 











Axis oF REvotuTION 


Gin @ 





Fic. 3. Nonlinear capacity, neck shaped vessel. 


In tanks of the form shown in Fig. 3, a is positive, and for the form shown in Fig. 4, 
a is negative. This expression also will give the approximate volume of many other 


tank forms. 
If a system consists of tanks, the volume of which may be expressed in the form of 


Eq. (2), the differential equation which defines the behavior of the system may be 


written as: 


C; ay. = iy + A2Y2 + G3 » (3a) 
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d 4 
C, dt (Yo + QY>) = by + b2Y2 + b; ’ (3b) 


where a, b, c with subscripts are all parameters. 

The purpose of this paper is to indicate a solution which describes the behavior of 
such a system. However to facilitate the analysis, it will be convenient to make certain 
transformations and to discuss the problem in terms of trajectories in the phase-plane. 
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Fia. 4. Nonlinear capacity, barrel shaped vessel. 





By solving for y, in Eq. (3b), differentiating and substituting the value of y, and 
dy,/dt into Eq. (3a), it can be shown that: 


2 


d d d 
mM, =3 (Yo + aty>) + mz = (y2 + ay) + Ms = Y2 + My. + m; = O, (4) 
dt dt dt 
where m’s are all parameters. 


By a suitable transformation, the differential equation, Eq. (4), may be reduced to 


the form: 


3 d , dY pes 
qe (¥ +BY) + Mag (¥ + BY) + Mig + Y =r. (5) 


The quantities 8, M, , M, , and r are all parameters. This equation differs from the 
standard nonlinear differential equations for which solutions are available. 
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If V is defined as: 
P adY 
ea 


then 
4 iy 4 BY’) = (14 38Y9V, 
do 


. i y 73, ? da ‘ 72 r 
qo? (Y + ,pY’) = 13 ay (1 + 36Y’)V]. 
Thus, from Eq. (5) 


dV ~ 


dV _ M, ) r— Y(1 + 68V") 
i aa 


i+ 38Y7, — (6) 





-(m, va+36Y) ’ 
where V and Y are the phase-plane variables. 

A discussion of the solutions of the key equation (6) in its general form will be 
temporarily postponed, and a special case where M, = M, = 0 will be considered, for 
in this case a solution in closed form can be obtained. 

2. An explicit solution for the special case, VM, = M, = 0. 

When M, = M, = 0, note that Eq. (6) takes the form, 


V dV rdY Y dY 


~J 








1+66V’ (1+ 36Y7\(1+ 68V’) 1+ 38Y? 

Equation (7) does not appear to be integrable because the first term on the right- 

hand side of the equation is a function of both V and Y. However, by introducing a 

new function ¢ and rearranging the equation, the explicit solution of Eq. (7) may be 
found. 








Define 
' r dY 
At} = i i (1 + 3Y7)(1 + 68V%) - 
Also, from Eq. (7), 
= " 12VdvV ¥ 128Y dY / 
oY) = . 1+ 6sv?* i 1+ 3gY*’ m 


where Y, and V, are the initial values of Y and V respectively. Equation (9) is readily 
integrable and 

rom Lt BBV (1 + 36Y)" 

oY) = nT 6avi (1 + 38¥2) - 

Let ¢, be the initial value of @. Using Eq. (10) the value of ¢, may be determined; in fact 


O(Y)acy-v.) = % = nl = 0; 


(V=Vo 


thus 





e-_at 6BV*)(1 + 3BY*)* 


Q , (11) 
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where 
Q’ = (1 + 68V5)(1 + 3BY%)’. 


Combining Eq. (8) and Eq. (11), the variable V may be eliminated and 





6 (1+ 38Y *) dy 
Y) = 12 ( 
Differentiating Eq. (12), 
db _ 159, € “(Ul + 38¥") 
qy = 126r ay (13) 


Now, the variables are separable and thus by integration 


¢ Y 
ef =a +ary] 
or 
o= 14 ZY + BY) — (Yo + BYD). (14) 


Combining Eqs. (11) and (14), the solution of Eq. (7) may be obtained: 
(1 + 68V*)(1 + 38Y’)’ = Q* + 126r[(Y + BY*) — (Y. + BYo)] 








or 
2 _ Sia (¥? + 1.58¥) + OY + BY?) 
where 
S’ = Yo + 1.58Y> + Voll + 36Y0)? — 2r(Yo + BYo) (16) 


is a function of the starting conditions. Equation (15) together with Eq. (16) yields the 
solution of the differential Eq. (7) in closed form. 

3. Bounded and unbounded solutions for r = 0 when M, = M, = 0. Bounded 
and unbounded solutions in the special case under consideration will be investigated. 
For any real value of Y and any positive value of 8 it may be seen that (1 + 38Y*)? > 1 
and from Eq. (11), setting r = 0, that 


V7? < S? —(Y’? + 1.58Y'), 


or 


V? + Y? + 1.58Y* < S’, (17) 


where 
S? = Yo + 1.58¥o + Vall + 38Y%)’. 
Thus V and Y are bounded in the phase plane as shown in Eq. (17) where the left- 
hand quantity has to be smaller than a positive real number S’. 


However, for a negative value of 8 the situation may be quite different. A plot of 
the solution of Eq. (7) on the phase-plane from Eq. (15) and Eq. (16) will demonstrate 
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this. It is clear that the solution of Eq. (7) depends on the starting conditions in the 
phase-plane. A numerical example where £ is negative will illustrate this feature. 


Assume 
——- = 0 
30° *™™ 
then 
Vy? = SS ~ 7 + 02" 
(1 — .04Y") 
These phase plane trajectories are plotted in Fig. 5. At S° = 12.5, the value of | V | is 
3.535 for every value of Y, except | Y | = 5.00. At | Y | = 5.00, V is indeterminate. 


The system is bounded provided the initial point (Y> , Vo) falls within the rectangular 
box | Y | < 5,| V | < 3.535. An initial point falling outside the box gives rise to an 
unbounded trajectory. 

For example, if the starting conditions are given as 


Y, = +4, Vo = +4 


we can see that the system is not bounded. 

4. The effect of changing the parameter 8. The choice of any particular value 
for 8, for example 8 = + 4/300, in no way limits the generality of the treatment, for 
the characteristic of the solution for any other value of 8 can be obtained by introducing 
a scale factor, i.e., a transformation parameter may be used to make the coefficient 6 
equal to + 4/300. More precisely, if this parameter is denoted by L, let Y = LZ, then 
Eq. (5) may be written as 


eae — ee a dZ ? r 
£2 + BL’Z) + M.— (2 + BD) + M, H+ Z=L, 
do do da L 
or 
f+) + M2 (24+) 4M, 2 4Z=R, (18) 
de doa do 
where 
t= ,6L* and R= =. (19) 


Hence ¢ may be made equal to either + 4/300 or —4/300 by a suitable scale factor 
L’ and R can be calculated correspondingly. The solution of Eq. (5) will be Y = LZ, 
where Z is readily obtainable by comparing Eq. (18) with Eq. (5). These remarks apply 
not only to the general equation, but also to the special equation arising when M, = 
M, = 0. 

5. Bounded and unbounded solutions for r ~ 0 when M, = M, = 0. In Sec. 3 it 
was shown that a solution would be bounded when r = 0 if the initial values (Y, , Vo) 
fell within a certain rectangle (see Fig. 5). When r ~ 0 this bounded region changes 
shape. Indeed the bounded region for Eq. (7) with a negative nonlinearity has a horse- 


shoe shape. 





| 
| 
| 
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Fic. 5. Negative nonlinearity for 8 = —4/300 andr = 0. 


To demonstrate this, Eq. (15) may be re-examined, keeping in mind that a bounded 
value of V is required as 1 + 38Y? approaches zero. In this case the numerator of the 
righthand side of Eq. (15) also approaches zero 


S? — (Y? + 1.58Y*) + 2r(Y + BY’) > 0. (20) 
This is the criterion for a boundary of the bounded region. 
By substituting Y* = —1/38 — e into Eq. (20), we have 
eee ae (-4)" 21 
‘ = 68 3 38 + €i ; ( ) 


where ¢ and e, are all infinitesimals. 
Since 8 is negative, B’ may be defined as 6’ = — 8. Using this value, 6’, and the 
value of S* from Eq. (16), Eq. (21) becomes 


2 _ (1/68") — (4r/3)(1/38")"* + 2r(Yo — B's) — Ye+ 1.58" (99 
‘ (1 — 36’Y,)° : 





The above equation determines the bounded region for Eq. (7). 

It is interesting to note that if Yi = 1/36’, Vo is indeterminate. However, the value 
of V? can be found as Y5 — 1/38’ by differentiating the numerator and denominator 
of Eq. (22) with respect to Yo ; in fact 


2r —_ 68'r YG - 3 Yo + 68’ Y> Yo Ts (23) 


wen Mi -WIO-@T) ~ wY. 
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If the solutions to Eqs. (22) and (23) are plotted for r = 0.1, 1.0, 2.0, 4.0 at 8 = 
— pf’ = —4/300, the boundary is seen to be of U or horseshoe shape as shown in Fig. 6. 


























-4 
Fig. 6. Bounded regions for different r, 8 = —4/300. 


At a particular value of r, the solution is bounded when any set of starting conditions 
lay within the horseshoe region; otherwise it is unbounded. The phase plane diagram 
for Eq. (15) is also given for r = 1.0; see Figs. 7 and 8. 


V 






, 
fo 











Fic. 7. Negative nonlinearity for 8 = —4/300 andr = 1.0. 
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Vv 
+ 











“6 


Fig. 8. Positive nonlinearity for 8 = + 4/300 and r = 1.0. 


6. The period and frequency when MV, = M, = 0. The period of a bounded solu- 
tion of Eq. (7) is in general an elliptic integral. This can be seen by rearranging Eq. 


(15) in the form 





.. ff .. ... 843 a7 
i) iin i) [S? — (¥? + 1.58Y) + 2r(¥ + BY] (24) 


If numerical values are assigned to the parameters, the integral of Eq. (24) may be 
evaluated by referring to any standard text on elliptic integrals. 
If r = 0, the general solution can be found by considering 


fv +ey)+¥=0. (25) 


Expressing o in terms of Y, Eq. (25) may be written in the form 


‘ . . ere 
I sales I [S?— Y? — 1.5eY‘]'?’ (26) 





where 
S? = ¥34+ 1.58Yo + Va(l + 38Y%)’. 
The solution of Eq. (26) is a combined elliptic integral of the first and second kind. If 


8 is negative, then the quarter period will be 


- S 2S 
p = (%2 - 28)xw + Sw, (27) 
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where 





-2 _ 1— [1 — 66’S’]'” E _ s" 
Y, =, ka || , 2 
38 Ss (28) 
Y,, being the maximum value of Y at V = 0, if bounded and p’ = —8. K denotes the 
complete elliptical integral of the first kind. E denotes the complete elliptical integral 
of the second kind. The general solution of Eq. (26) will be in the form 


‘ 28 2S , 
o= (<= -_ 280, k) + 7 E(6, k), (29) 
where 
@ = sin’ J . (30) 
re 
If 8 is positive, the quarter period will be 
S 2 1/2 
p= +[ (8) -- | [2E(k) — E<k)], (31) 
where 
elt [i+ sy? L E =," 
Fe . 38 ? k 7 9S? ie, ee ’ (32) 
Y,, being the maximum value of Y at V = 0. The general solution of Eq. (26) will be 


in the form 


s\? 1/2 Ix : lx 
¢= | 2-8 ) - | | 2n(% k) — 2E(6, k) — K(E 2 k) + F(86, 0], (33) 


where 


a 


6 = cos" 7 (34) 


The frequency may be calculated for both negative and positive 8 as 
a/2 
a 
p 
The frequencies for various values of Y,, and for 8 = + 4/300 are plotted in Fig. 9. 
It is interesting to compare Eq. (25) with the Duffing’s equation, 


2 


a. X+X+,)X*=0, (35) 
where the period is an elliptic integral of the first kind only. If 8 and \ are small, the 
approximation \ = — 8 = #’ may be used. 

Therefore the negative nonlinearity of Eq. (25) corresponds to the positive non- 
linearity of the Duffing’s equation, and vice versa. This can be seen as the positive 
branch curves towards the lower frequencies just as the negative branch does for the 
Duffing’s equation. However, the frequency will be entirely different as starting ampli- 


tude gets higher. 
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7 ? 
| 7 (Y + BY’) + ¥ =0, 
"7 Y,, = Max. Y, 
| _ 7/2 
12 aaa P ’ 
P = Quarter Period. 
5£ 10 
wW 
a 
> 
=8 
= | 
a } 
61 





‘5 6 7 ‘8 ‘9 1-0 re 
FREQUENCY W 


Fic. 9. Frequency amplitude response. c denotes cut-off point at ¥m = 5.0. 


For a very high amplitude the frequency slows down for a positive 8, while there 
is no bounded solution for a negative \ with a very large amplitude in the Duffing’s 
equation. For a negative 8, the amplitude is restricted and not more than (1/36’)'”, 
while there is always a frequency for any positive \ at high amplitude for Eq. (35). 

7. Solution with damping terms. With all these details in mind for the character- 
istic of solutions of Eq. (7) in the special case 7, = 1, = 0, let us turn now to the 
characteristic of solutions of the general equation (6). Note first of all that the slopes 
of trajectories in the phase plane, as specified by Eq. (7), must be decreased by an 
amount A/, + [M,/(1 + 38Y°*)] if the slope of the trajectories is to be that specified by 
Eq. (6). In other words, when M, and M, are not zero, damping occurs and the slopes 
in the phase plane must be adjusted accordingly. Conceivably, then, one would begin 
with a solution of Eq. (7), sketch in a number of isoclines (i.e., a number of tangent 
line elements); adjust the slopes and draw new isoclines for Eq. (6), and then sketch 
in the trajectories for Eq. (6) and have a picture of the performance in the large for the 
general case. 

For the general case, however, it was more convenient to use an analog computer 
which sketched the trajectories for the following cases: 
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B - 300 ’ 

M, = 0.2,M, =0 in Fig. 10 

= 0 in Fig. 11 


r = 2.0 


— = 

~ ~ 
~ ~ 
Il Il 

S — 
~ . 

~ 

i) 
~ 

— — 

= ~ 
2 > 
ll | 


0.2 in Fig. 12 


M,=0, M, = 1.0in Fig. 13 


Note that if 6 > 0, the quantity 7, + [M,/(1 + 38Y°)] representing the change in 
slope is always positive and the effect will be the creation of inward spirals toward the 
value r. If 8 < 0, the change in slope is also positive for bounded solutions since the 
quantity (1 + 38Y”*) is always positive for bounded solutions without damping. More- 
over, (1 + 38Y’) <'0 is not possible in a physical situation. 











Fra. 10. Computer solution, nonlinear capacity, for 8 = —4/300, r = 2.0, Mi = 0.2, and M, = 0. 


8. Bounded solution when M,, ~ 0 and M, = 0. It can be shown from Figs. 12 and 
13 that the higher the value of M,, the larger the bounded region. This bounded region 
for any particular value of M, may be obtained by numerical integration of Eq. (6) 
provided that the initial values of Y and V are known. It is logical to consider the upper 
right corner of the bounded region in Fig. 6 as a point through the integral curve. For 
any given value of r this point is determined by Eq. (23) as shown in the following 
example, where r = 2, Y. = 5.0, and V. = 2.7387. 
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Fic. 11. Computer solution, nonlinear capacity, for 8 = —4/300, r = 2.0, M,; = 1.0 and M, = 0. 
V 








Fic. 12, Computer solution, nonlinear capacity, for 8 = —4/300, r = 2.0, M, = O and M, = 0.2. 
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‘ae 











Fic. 13. Computer solution, nonlinear capacity, for 8 = —4/300, r = 2.0, 4, = Oand M, = 1.0. 


A closer examination of Fig. 6 and Eq. (23) indicates that this corner point obviously 
has two distinct slopes—one being infinite, the other being finite, depending upon the 
direction of approach to the point. To determine the finite slope Eq. (6) is rewritten in 
the following form: 

; dV r— Y(1 + 6BV") 
lim =s=—-M,+ lm —;--ScR 
Y ‘ rszay20 ©606-: V(1 ~ +s 3B Y°) 


1+38Y2-—0 d 1+3 


By denoting the limit values of V and Y as V, and Y, respectively, it can be shown that 


1V 72 ee ; 
(¥) = —(1 + 68V(186V.Y.)' — M,/3. (36) 
The limiting value of the second derivative of V with respect to Y can also be expressed as 


_(@V\ _ {ey deg , (dV : me 7 
(<4) = V? ay). + (1.25Y>' + 0.50M,} (FF) + 0.25M,Y;'. (37) 


For the same example with r = 2 and M, = 1.0, one obtains 
I ’ 


\ 


dV ’ : (gv ox 
(40) = —,2116, eI = 0251. 
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To start the integral curve the Taylor’s expansion with three terms is adequate. 
incidii d v) h? (¢ ‘) 
an 7 } naa 4 — | —- te eee 
Ata (%y a 2 \dY’/. 


where h is the increment of Y. The mth approximation of f may be defined as 


— Y,(1 + 68V2) 
V.(1 +38¥2) (38) 


dV . r 
—_—_ = = of + 
dY fal ) ee 





The above formulae together with the following Simpson’s rule constitutes the 
method of successive approximation in numerical integration. 


V.(a + 2h) = V(a) + : (f(a) + 4f(a +h) + fala + 2h). 


TABLE 1 

Bounded region with damping term M, forr = 2, M, = 1.0, Mi = 0 
r V dV/dY y V dV/dY 
5.0 2.7387 — .2116 —5.0 — 4.1833 — .4130 
4.8 2.7815 — .2166 —4.8 —4.2671 — .4254 
4.6 2.8253 — .2216 —4.6 —4.3535 — .4377 
4.4 2.8702 — .2283 —4.4 —4.4425 — .4542 
4.2 2.9165 — .2336 —4.2 — 4.5350 — .4680 
4.0 2.9637 — .2401 —4.0 — 4.630 — .484 
3.6 3.063 — .253 —3.6 —4.830 —.519 
3.2 3.167 — .269 —3.2 —5.045 — .558 
2.8 3.278 —* 286 —2.8 —5.277 — .603 
2.4 3.396 — .306 —2.4 —5.528 — .655 
1.6 3.662 — .356 —1.6 —6.102 — .787 
0.8 3.972 — .429 —0.8 —6.802 — .976 
0 4.359 — .541 0 —7.689 —1.260 
—0.8 4.857 — .727 0.8 —8.865 —1.721 
—1.6 5.565 —1.073 1.6  —10.540 —2.545 

—2.4 6.677 —1.809 





The same technique can be applied to the lower left corner of the bounded region 
in Fig. 13. The results are listed in Table 1 and plotted in Fig. 14. This is to say that , 
two integral curves determine this bounded region. 

On the other hand, for M, = 0.2 and M, = 0 there exists only one integral curve. 
It is sufficient to determine the bounded region by this simple curve alone. 

9. The general bounded solutions. The value of M, (1 + 38Y*)~’ approaches 
infinity as (1 + 38Y*) becomes zero. However, the corner point with two distinct slopes 
shifts to a new position by letting 


—Y(1 + 68V2) — M,V. +r =0, 


from which is obtained 


M, ( M, J r— vr." 
= = ek Ta a 
Ve = —Toay, * | iosy,) + eay. | ° (39) 
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Fig. 14. General bounded regions for 8 = —4/300 andr = 2.0. 


The limiting value of the finite slope and second derivatives are also obtained as 





(7) = —(1 + 66V? + 68V.Y.M,)(188V-Y, + M,)"', (40) 
*) - ey ee = , “g 
-( - |v: (Sy) + (.25¥2" + 0.50M,~) a) + 0.25M.Y: si 


[1 + M,/248Y.V.J". 


The same method is used in numerical integration. The results are plotted in Fig. 
14 for r = 2, M, = 1.0, and M, = 0. Figure 14 indicates that the bounded region with 
M, = 1.0 and M, = 0 is bigger than that with M, = 1.0 and M, = 0. Thus it can be 
concluded that the M, term yields more weight in damping than the M, term. For 
the same type of damping, the lower the value of the M, or the M, term, the smaller 
the bounded region. This can be shown by comparing Figs. 10 and 11, or Figs. 12 and 13. 

For M, = 0 and M,, being finite, there are always some unbounded solutions of the 
differential equation. This is also true for M, = 0 and M, being finite, except the value 


of V. being higher. 
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10. Further physical applications. An electric circuit with linear resistance and 
capacitance in series with a coil having inserted magnetic material can serve as a good 
example for illustrating more explicitly the results in physical terms. The coil has N 


























(A) 
B t_ SATURATION 
H 
(b) 
Ee: E. 
C4 Rt q 
: i 
(C) 
E 








nr 





tease © 
Fig. 15. Nonlinear circuit: a) saturation, b) circuit elements, c) E denotes V(1 + 38Y?) as 1 + 38Y2-—0 


turns with a magnetic flux ¢, thus the emf across the terminals of the coil as shown in 
Fig. 15(a) is 
d 
& = dt (N@). 


The differential equation can be written as 


dt 


To solve the above equation it is necessary to determine the relation between 
@ and q. If the flux density is B and the magnetising force is H, the relation between B 
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and H before saturation may be approximated by 
B = K(H — aH’), (43) 
where 
K = constant, a = nonlinear parameter. 
Let u be the incremental permeability of the material, then by differentiating Eq. 
(43) 


dB 


deal aa K(1 — 3aH’). 


If H is small, u is approximately equal to K as its limiting value. The flux density 
is saturated at 


dB 


dH = wl — 3aH”’) = (). 


Thus the B — H relation can be defined only within the region 
1 — 3aH’ > 0, 


or 
H® < 1/8a. (44) 


Substituting the value of H = Nil™' and B = ¢A™. 


¢=S'Ni-— yi), (45) 
where A is the cross section area of the material, 1 the length of the magnetic circuit, 
i the current, S = 1K~'A~* the incremental reluctance, and y = aN’1~? a nonlinear 
parameter. 

Thus, 
a P , 
& = dt (N¢) - Lat had y), (46) 


where L = S"'‘N’ = inductance for the linear portion. 
By substituting e, in Eq. (42), the differential equation becomes 


* ne 3 : . 
L5,@ — vi) + Ri +! = B,. (47) 


The differential of the above equation will yield, 


ee di,i_ dK, 
LG - 1) +RE+- = =. (48) 


For a step change of E, dE,/dt = 0 it can be seen that this is reduced to 


2 


4. (@¢ — yi’) + (RL“'?C'”) a +i=0, (49) 
do da 
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where 
o = (LC), 
or 
2 
fy+eY)+M, 2+ Y¥ =0, (50) 
do dco 
where i = JY, 8 = —yI’* = —4/300, M, = RL™'?C™’, and J is a reference current. 
For a control device the system sometimes follows a ramp type signal defined by 
E, = tt, (51) 


where 7 is a constant. 
The system will have a differential equation 


ae “tf oe 

£66 -4 4+ GORE ties (52) 
do dt 

with the starting conditions at t = 0,7 = 7% , di/dt = 0. Equation (52) can be expressed 

in nondimensional form as 


-— rs dY 
qe (¥ + bY) + Moe + VY =r, (53) 





where r = rl‘ ato = 0, Y = Y, = ipl", dY/do = 0. 

By knowing the values of L, R, C, y, r and 7 it is possible to compute the parameters 
L, M,.,rand Y, with the fixed value of 8 = —4/300. To determine whether the solution 
is bounded it is a simple matter to examine the starting conditions in the phase plane 
curves for the above nondimensional parameters. 

The definition of boundedness as given in previous discussions is entirely mathe- 
matical. The unbounded solution was interpreted as the value of V not bounded. How- 
ever, close examination of Eq. (15) indicates that the product of V and (1 + 38Y7’) is 
bounded. This can be written as 

V(l + 3BY’) = [S’ — (Y? + 1.58Y*) + 2r(Y + BY*)]'”. (54) 

The curves in Fig. (15) give the values of V(1 + 38Y’) as (1 + 38Y’) approaches 

zero. It can be proved that the damping terms M, and M,, do not affect these values. 


Equation (6) may be rearranged as 





ies 5 ie dl ond r— Yl + 68%) . 
de ‘¥ + BY) = Vil + 38Y) = Gray tM. + Ma+eery = 


As (1 + 38Y”’) approaches zero, V* and dV /dY are of higher order than (1 + 38Y’)~". 
The importance of introducing the above analysis is that the voltage e, is always 
finite at the saturation flux. This is because e, can be expressed in the form of Eq. (55). 
Physically the rate of change of Y, (i.e.,V), rises very rapidly as Y approaches (1 —1/38)'”, 
but the value of Y is only defined within this region. 
To carry the physical problem to a new region, the final condition of the nonlinear 


differential equation is of interest. 


s = at, — Ri, —¢,, (56) 
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where 7, and e, are final values of 7 and e, respectively, and ¢, is the time duration for the 
current 7 becoming 7, . The new region has a horizontal line on the B vs. H curve after 
the flux is saturated, ¢ being constant. It is therefore concluded that the value of e, 
suddenly drops to zero outside the nonlinear region. The system will then become a 
linear RC circuit with the same ramp input 


dg, g_ 
Ra a = gt, (57) 


The starting conditions of this linear differential equation are the same as the final 
conditions of the nonlinear differential equation. 

11. Conclusions. The solution of the system of differential equations (3a) and (3b) 
can be conveniently represented in the Y V-phase plane and various sets of graphs made 
to show the effect of varying the parameters 6, r, M, and M,. For a positive 8 the system 
is always bounded and the frequency tends to be slower than a corresponding linear 
case. For a negative 8 with no damping the system is bounded within a horseshoe region 
and the frequency tends to be faster. If the system is bounded without damping, it is 
always bounded with damping. The effect will be a spiraling towards the value r in the 
phase plane. 

In a physical system an unbounded solution may be re-defined by extending to a 
new region beyond its boundary. The analysis can be applied to liquid level controlled 
systems or circuits with flux saturation. 
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EXPANSIONS OF THE IRREGULAR COULOMB FUNCTION* 


BY 
FRANK 8S. HAM** 
Department of Physics, Harvard University 


Abstract. Convergent and asymptotic series in the energy are derived for an irregular 
Coulomb function appropriate to both attractive and repulsive potentials. The co- 
efficients in these series are Bessel functions. For repulsive potentials the derivation 
yields series in agreement with those given by Breit and others. Some of the series for 
attractive potentials have been given by Kuhn. The present general derivation and 
some of the series are, however, new. Tables of the coefficients in the series for the 
attractive potentials for values of the parameters of interest in solid state energy band 
calculations have been prepared and are available. 

I. Introduction. Breit and Hull’ and Abramowitz’ have shown that the irregular 
Coulomb function G, can be expanded in an asymptotic power series in the energy which 
is valid for repulsive Coulomb fields. Jackson and Blatt® have given a convergent series of 
similar form for G, , and their derivation can evidently be extended to higher L. In the 
present note we shall derive convergent and asymptotic series for an irregular function 
appropriate to both attractive and repulsive potentials and shall obtain the results of the 
above authors as a special case. A review of the numerical treatment of Coulomb func- 
tions has been given recently by Fréberg (Rev. Mod. Phys. 27, 399 (1955)). 

The investigation leading to these general results was undertaken in order to establish 
formulas for the convenient numerical calculation of an irregular function for attractive 
Coulomb fields and positive and negative energies of small absolute value. Such calcu- 
lations are an essential preliminary to the theoretical investigation with the “Quantum 
Defect Method’’* of the energy band structure of the alkali metals, a program of research, 
which is being continued by Brooks and the author. Therefore, we shall use a notation 
in this paper which is most appropriate to work with an attractive Coulomb field, but 
in an appendix we shall show the connections with the notation used by Breit and 
others in their work with a repulsive field and shall show that our results agree with 
theirs. 

Kuhn’* has already given the asymptotic series that will be derived in this paper, 
but his presentation is marred by the fact that he did not realize that these irregular 
series do not converge. Our convergent series are, however, new. We shall assume that 
the reader is familiar with Kuhn’s paper. 
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A set of tables of the coefficients in the series for the regular and irregular functions 
for the attractive Coulomb field has been prepared for values of the parameters of 
interest in the solid state calculations. These extend and correct the tables given by 
Kuhn’ and are accompanied by a more detailed discussion of the formulas than can be 
given here. They contain also tab!es of various quantities which appear in the formulas 
and which will not be written out explicitly here. These tables are now available and 
may be obtained from the author or from the Division of Applied Science, Harvard 
University’. 

II. Convergent series representations of Coulomb functions for non-integral values 
of (2L + 1). The hydrogenic radial wave equation in Rydberg units for an attractive 
potential takes the form’ 


"U 1 
7 +{-4+2- 
rr n 


= Ibo 


L(L + Dp U =0, 1) 


r° 


where r is the radial coordinate, U the radial part of the wave function multiplied by 
r, L the angular momentum quantum number, and E = —1/n’ the energy. Yost, 
Wheeler and Breit’ and Kuhn’* have shown that solutions of (1) exist which have the 
form 


US” (2) _ Yn * Ui), (2) 
k=0 
where z = (8r)’””. If we set 
LP@ = @/2)Vi"@), (3) 
we find that the V;”’(z) satisfy the simultaneous differential equations 
V7. Vie 2) = 0, (4a) 
ViVi @ = @/2)'Vii®, (4b) 
where 
a 
vi=e(S :) + (4 L) at — (OL + 1)°. (5) 


Kuhn has shown from the recurrence relations satisfied by the Bessel functions that 
if C,,(z) is a linear combination of the Bessel function J,,(z) and the Weber function 
Y,,(z) 





C,,(z) = AJ,(z) + BY,(z), (6) 
A and B being independent of m, then 
2L + 2 + T4Tq /2)** - l ¢ a+3 
v2 A(2 + 9) (2/2 "Corss+e(2) — 434+ 4) (2/2) Carsendl} (7) 
= (2/2)**Corsisel2), 


and solutions of (4) may be generated if we choose 
Vo) = Cozsil2). (8) 


*F. S. Ham, Tables jor the calculation of Coulomb wave functions, Technical Report No. 204, Cruft 


Laboratory, Harvard University 
7F. L. Yost, J. A. Wheeler and G. Breit, Phys. Rev. 49, 174 (1936) 
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In particular, if in (6) we choose A = 1, B = 0, and if we then determine V{”’(z) from 


(7) for k > 1, so that no V;”’(z) for k > 1 when expanded in powers of (z/2) contains 


the power (z/2)’”*', we obtain the series 
2 3k 
n \ / 0 (Lin —2k (L) ‘ 
Jizai(2) = (2/2e) US) = Sn™ YS a (e/2)* Jarsree2)- (9) 
k=0 q=2k 


It is shown in Appendix I that this series converges absolutely and uniformly for | z| < R 
and || > mn , where R and mn, are arbitrary positive numbers. 
Similarly, we may show that 


—2L + 1(2)""¢ a 1 (3)"" 
742 + q) 9 C_ (2041) +¢+2(2) 4(q 3 9 C90 +1) +e+2(8) (10) 


> ers 
= (2) C2041) +e2)- 


Choosing Vj” (z) = J-.2n41)(z), we obtain the series 
«© Bk 
n —2 (L) 
J (au +1)(2) = dn , } by e(2/2)°S — (2041) +0(2)> (11) 
k=0 q=2k 


in which the }{”) are determined from (10) in the same manner as the a‘? are obtained 
from (7), so that (z/2)~°**” occurs only in the term V4”? (z). In the following discussion 
we shall always mean by a‘’) and b>, the coefficients appearing in these series, de- 
termined as described above. The proof that (11) converges absolutely and uniformly 
for 2z| < Rand!n| > np is carried out in a manner identical with the proof of the 
convergence of (9). 

Since both (2/2)J3,.,(z) and the regular Whittaker function® M,,,,;(2°/4n) satisfy 
(1), it follows by comparison of the coefficients of (z/2)°*? in the series expansions of 
both functions that° 


L+1 


Jets) = | ey |Me-esrale/An) (12) 
Similarly, if (2L + 1) is not an integer, we may show that (2/2)J* .21+1)(z) is related 
to M,,,-1-,(2°/4n) by a formula obtained from (12) by replacing L by —L—1. Hence 
if (2L + 1) is not an integer, (9) and (11) are two independent solutions of (1), both of 
which are satisfactorily convergent in both z and 1/n. Kuhn has proved the existence 
of two such series solutions’ and has exhibited (9); however, instead of (11) he supposed 
that the series obtained from (9) by replacing J,,(z) by Y,,(z) converged, whereas in 
fact the resulting series diverges, as we shall show. 
If (2M + 1) is a positive integer, the function M,,_y_;(2°/4n) is not defined, and 


Q 
we may show that 





+/ _(9 
ties Shoes » SOS Dee SS pe, an. (13) 
LM n°” T(n — M) 


8E. T. Whittaker and G. N. Watson, A course of modern analysis, 4th ed., Cimbridge Univ. Press, 
1952, p. 337 


’G. H. Wannier, Phys. Rev. 64, 358 (1943) 
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Hence in order to obtain a second solution of (1) which is independent of J2,,,(z) for 
all values of 2L + 1, we define’ 


Nivel) = Ez a =a 2 Fiss(@) cos (2L + 1)" — Ja: © |(ataas) (14) 
and understand that, for integral values of (2Z + 1), N2z.,(z) is defined by the limit 
of the expression on the right of (14) as (2L + 1) — integer. 

III. Series representations of N2,,,(z). To derive a convergent representation of 
N21+1(2) in @ power series in (1/n”) we may substitute (9) and (11) into (14), since 
both of these series converge absolutely for all L. However, for non-integral values of 
(2L + 1) we cannot expand the coefficient of J2,,,(z) in (14) in a convergent power 
series in (1/n”) because !'(n + L + 1) has a pole whenever (n + L + 1) is a negative 
integer or zero. We therefore write (14) as 


Tm+L+1)_ 1]. 


N2,.,(2) sin (2L + 1)x = | 
rF n®“*'T(n — L) 


| rs > Q¥.e(2/2)"Jarsr+e(2) cos (2L + 1)x 


ne (15) 
2 Sk 
+> nn a\?(—)"(z/2)"Jors14¢(2) cos (2L + 1+ g)x 
k=0 q=2k 


3k 
(L) ¢ 
~ z. by. e(2/2)°S ~ 12241) a}. 


q=2k 


The k = 0 term of the second series in (15) is seen to be sin (2L + 1) Y2,,;(z), since 
ais? = bi? = 1. Since the first series in (15) is simply a multiple of J2,,,(z), and because 
both (z/2)J2,.,(z) and (z/2)N2,.,(z) satisfy (1), we see from the form of (15), Eq. (7), 
and the above value of the k = 0 term that the k = 1 term of the second series in (15) 


is determined by (4) to be 
n *la\2(z /2)*¥ “et+3(2) + .3(2/ 2)° Yo1+4(2)] sin (2L + l)x (16) 


plus some linear combination of the solutions of the homogeneous equation (4a), J2.,(z) 
and J_,,-,(z). We can determine this latter combination by using the Bessel function 
recurrence relations to re-express the terms involving (1/n’)J_2,-1,,(2) in the second 


series. We see that we shall obtain 
— b5°2(2/2)?J -on-142(z) — 0 9°3(2/2)° J -22-143(2) 
= —@'‘)(2/2)*J_or-3(2) + a'\9}(2/2)* J -21-a(2) (17) 
A(L)J-21-:(2) — BUD) J214:@) 


We need not carry this out explicitly in order to obtain the polynomials A(L) and B(L), 
however. Because the recurrence relations relate Bessel functions of orders differing by 
integers, if 2(2L + 1) is not an integer’® we shall not obtain J2,,,(z) in an expression 
for J_21-,.,(z), so that B(L) = 0. We may obtain A(L) if we notice that when 2L + 1 
equals an integer, 2M +- 1, the left side of (15) vanishes for all n and that 


10We make use of non-integral values of 2(2L + 1) in order to prove B(L) = 0. Both sides of (17) 


being continuous functions of L, the conclusion holds for 2(2L + 1) = integer as well. 
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(M+1] 
nM +n +1) _ BAM) (18) 
n IT'(n — M) po On” 


the b,(M) being polynomials the calculation of which is outlined in Appendix II, and 
[M + 1] denoting the largest integer less than M + 1 (i.e. [M + 1] equals M if M is an 
integer, (M + 4) if M is half an odd integer). Therefore the coefficients of (1/n”) on 
the right of (15) must vanish when L = M, and since J_2y-,(2) = (—1)?”*"Jams.(2), 
we find on combining all terms in 1/n’ in (15) when L = M and using (16), (17), and 
(18), that A(L) must differ from b,(L) by at most a function of L that vanishes for all 
M. Since A(L) is itself a polynomial, however, we have that A(L) = b,(L). 

Continuing this argument to higher powers of (1/n’), we find that the desired form 
of the series is given us immediately: the term in n~™ is determined from that for n~?“~” 
by the use of the differential equations (4) and (7) except for a multiple of n™™*J_21-,(z); 
this multiple is determined from the requirement that this term cancels for all integral 
values of 2M + 1 the term —n~™*b,(M)Jo.,(z) in the first series of (15). 

We now wish to use this information to rewrite (15) in a form such that when L = M 
each term of the series vanishes. We can’t do this, once and for all, for arbitrary M because, 
as we shall see, we should have to make use of a divergent asymptotic series which we 
want to avoid. But if we are interested in a particular M we can conveniently proceed 
as follows. We know that the series (9) and (11) converge absolutely, so that we can 
add to the second series in (15) a finite polynomial in powers of (1/n’) multiplied by 
(9) and rearrange the resulting combined series in powers of (1/n”) without disturbing 
convergence. Choosing this polynomial to be >>!™\""' b,(L)n~”’, adding it in as described 
and subtracting it out in combination with the first series, we have finally on making 
use of the relation b)(L) = 1 


- {M+1] 9 a 
Nivel) = f (L+n+i)_ "ss bf Ot + De 5: 40 


nit 1TD(n — L) a0 «on™ J sin (2L + 1)x 





[M+1] 3(k—p) 


+ ps ae > b,(L) oH Ox) (2/2) *Vor+1+0(2) 
k=0 p=0 q@=2(k—p) (19) 


{M+1] 3(k-p) 


>< p> b(L) 2D) asny,q(2/2)*Yarerse(2) 
k= (M+1) 41 


p=0 a=2(k—p) 








7 > | b,(L) | i as!) —)"6/2)'T arson}. 


p-twaij+1 LSIN (2L + 1)r] 20-») 
We see that the series is now in a particularly convenient form for taking the limit 
L — M, where, we recall, 2M + 1 is a particular integer. Thus using (18) we can write 
for L— M 








_. cos(2L + Ir {r(L+n+1) _ bi) 
Gu(M,n) = lim sin (2L + 1)m ee a — 


1fd/rma+L+4+)D i _, d \ 
ie: soe te ee eee ha > ° b(L 
9 (f (Ets ite ») p=0 ¥ dL ( of )) 


aT 


LeoM (20) 
a 2 {re +2+) ym 4141) +¥m—- LD) -2n@)) 


n’’* Tin a L) 
{Af+1) 1 d 
-> >> OL 0) 


p=0 








36 FRANK 8S. HAM Vol. XV, No. i 


where ¥(x) = (d/dx) In I'(2). Also it is in general true that b,(M/) = Oif p> [M + 1]. 
Thus all terms in (19) vary continuously as L — M. 

We notice in (19) that (z/2)N2,,,(z) has been expressed as a multiple Gy(L, n) of 
the regular function (z/2)J3,.,(z) plus an absolutely convergent series in (1/n*) which 
we shall call (z/2)Q5,.,(z, M) and which is clearly an irregular solution of (1). Both 
(z/2)Nor4:(z) and (z/2)Q3,.,(z, M) satisfy a Wronskian relation with (z/2)J2,,,(z): 


[(2/2) Jo1+1(2)[(2/2)(d/dz)(z/2)N314,(2)] (21) 


— [(z/2)N21+:(2) [(e/2)(d/dz) (2/2) J314:(2)] = 2°/4e. 


i\< 


The relation involving Q2,.,(z, M) is identical with (21) if N is replaced by Q. A number 
of other formulas involving N2,.,(z) are given in Appendix III. 

An asymptotic series for N3,,,(z) in powers of (1/n*) for arbitrary L may be obtained 
from the convergent representation (19) if we make use of the asymptotic expansion 
(A3) (Appendix II), valid for | arg (n) | <  —A < 7, and rearrange (19), making use 
of the fact that we may multiply and add asymptotic series. Combining terms in 


Jozstee(2) and J_s; .(z) by means of the definition '’ of the Weber function, we obtain 
x x 3k 
= —2p — 2k L) 9) , ) 
Nir) ~ Do 2n-”b,(L) Don™ >, ab(e/2)*Vere1+e(2)- (22) 
p=0 k=0 q=2k 


Multiplying by the reciprocal expansion of (A3), we have also 


T'(n — L)n?**' oe ' : a ae “L) a f 
a tL ap sen@) ~ FF Oe VN iscecl- (23) 
4T In|— k=0 


q@=2k 


Except when (2L + 1) equals half an odd integer, both (22) and (23) may be proved 
to be divergent series® which are asymptotic representations in n, for | arg (n) | < 
x—A < 7, of the functions on the left of the equations. We note that (23) has the same 
form as the convergent series (9) but with Y,,(z) replacing J,,(z)."* 

IV. Acknowledgments. The author would like to thank Professor Harvey Brooks 
for discussing this work with him and for making several helpful criticisms. 

Appendix I. Proof of the convergence of (9). It follows from (7) that the numbers 
a‘) have an upper bound 6‘”’. From the series definition of the Bessel function J,,(z) 
we may show that 

J n(z) | < (| 2 |/2)"(m!) ; exp (| z */4). (Al) 


Since 2k < q < 3k in (9), by writing (9) in the form (2) we may show that, if /, is suf- 
ficiently large, then for all k > ky , 


"UL (e)n™™* | < B'” exp (| z |?/4)(| 2 |/2)?*****? | n |-™ [(2L + 2k!) '. (A2) 


If we replace | z| by R and |! by n, on the right of (A2), where R and ny are arbitrary 


positive numbers, then the series in k formed from the resulting quantities can be shown 
uWhittaker and Watson, op. cit. p. 370 
“The series (23) is one of those given by Kuhn in Footnote 5. We have also verified for L = 0, 
1, 2, 3 that Kuhn’s alternate generating procedure with Co(z) = Yo(z) and Ci(z) = Yi(z) yields the 
asymptotic series (22) for (z/2)N3,.,(z) for these values of L, for which b,(L) = 0 when p > L +1. 
Hence Kuhn’s tables for the irregular function are convenient if the asymptotic representation is suf- 


ficiently accurate. 
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to converge absolutely by the ratio test. It therefore follows from (A2), the comparison 
test, and the Weierstrass M-test that the original series (9) converges absolutely and 
uniformly for |z| < R,|n| > m. 

Appendix II. The polynomials b,(//). From the exponential form of Stirling’s 
series'® for I(x), we may show that for an arbitrary L the following asymptotic relation 
holds for | arg (n) | < r—-A < rwas|n|— @: 


Ta+Ll+)) ei > b,(L) (A3) 


2L+1_ 2p 
n?”’*'Tin — L) p=o 


the b,(L) being polynomials in L. In order to determine these, we may let L equal an 
integer M. Then from the relation T(x + 1) = xI'(x) we may evaluate the left side of 


(A3), obtaining 


M 
n> "(n? — M’)(n? — (M — 1)’) -:: (rn? — 1) = : b,(M)n-*’. (A4) 
p=0 
Clearly b,(L) vanishes if L equals an integer less than p. It also vanishes if L is a half 
integer less than p — 1. We see from (A4) that b,(/) for integral M is equal to (—1)? 
times the sum of all different products of the squares of p integers less than or equal to 
M, no two integers in any product being the same. This may be written as 


M kp-1 kp-1-1 ka-1 
b(M=(-? > BD BB, DS #B.--- DA. (A5) 
kp=p kp—.=p-1 kp-a=p-2 kim 
Therefore 
M 
b,(M) = — >> k2b,_,(k, — 1). (A6) 
kp=P 


The relation (A6), the fact that b,(4/) = 1, and the formula 
M 
>> k(k + 1k + 2) «++ (k+7) = M(M + 1)(M 4:2) ---(M4+n+ (n+ 2) (AZ) 


may be used to calculate the b,(/) as polynomials in. M. Since the resulting polynomials 
are correct for all integral values of M > p, we can replace M by L in order to obtain 
the b,(L) appropriate to (A3) for general values of L. These polynomials are given in 
the reference of Footnote 6 for p < 7. 

Appendix III. Miscellaneous formulas. The following formulas are useful in manipu- 
lating the Coulomb functions. They may be derived from the definitions of the various 
functions and formulas given in Chapter 16 of Whittaker and Watson*. We should note 
that (A8) corrects an error in the relation of Example 2, p. 346 of that text. Also, (A9) 
was given incorrectly by Wannier® and corrected by Kuhn’. 


_ - 
aa oe! ) 0 Wn n+ 1s2(se'") 
. : (A8) 
TOL + 2) iem-t-vyy 
rin + L +1)‘ Wa, n+rsals)- 


M,, L+1/2\8) = 


((2L + 1) # negative integer]. 


BWhittaker and Watson, op. cit., p. 253 
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Wa1+1/2(2"/4n) = [T'(n a L + 1)n7*-* 2/2) J¢14+1(2) cos (n — L — 1) (A9) 
+ T(n — L)n(2/2)Niz4,(2) sin(n -— L—1)x]. 


@/2)Ninas(@) = (n*e'**)(1 a)" sin x(n — LEE + 1 — 0) Wor srra@/4n) (4 19) 
— T(L+ 147) cosa(n — L)W_,. 241;2((2"/4n)e**)]. 








; r(2L + 2) = mL oe" is, 
Moca) ~ en fs pony, All 
sum ~ Fit+i-0° “ trken * * (All) 
[((2L + 1) # negative integer; —r < arg (s) < 0.] 
(2/2)Noz.i(2) ~ (n-“e'"*)(1/n) 


[e**"*-' sin x(n — L)T(L + 1 — n) exp (—2?/8n)(z?/4n)" (A12) 
— (LL + 1+) cos x(n — L) exp (2’/8n)e**"(z’/4n)"]. 
[—r < arg (2*/4n) < O]. 


Appendix IV. Derivation of formulas for the repulsive field. Equations (9), (19), 
and (23) are valid for arbitrary complex values of the parameters L, z, and n. From (1) 
we see for attractive fields that if we take r and z positive real, then real (and for con- 
venience, positive) values of n correspond to negative energies, pure imaginary values 
of n to positive energies. Again from (1), we see that for repulsive fields and positive 
energies we may conveniently choose r to be real and negative and n and z to be negative 
pure imaginary numbers. We shall restrict our attention to values of Z that are positive 
integers or zero. 

The regular and irregular functions, F, and G, respectively, used by Breit and 
others'~*, have the asymptotic forms for fixed energy and large radius 

F, ~ sin [p — (rL/2) — 9 In2p + arg T(L + 1 + 1y)], 
Se (A13) 
Gi ~ cos [p — (rL/2) — n In2p + arg T(L + 1 + @y)]. 


Here in Breit’s notation p = | r|/|n|,7 =|n!,and2 =|z| = (8pn)'”. Comparing 
these with the asymptotic expansions of our functions for large | r|, got from (A11) and 
(A12), we find that 





T(L ‘i. 1 + in) Lg teste (een (z 
Qin’ ** \2 


F, ) zit (ia) (A14) 


and that 





(= weit. (—ia) = |TL+1+ 4m) | oe (1/m)fe" GF 1) — eG]. (A15) 


2 
From this latter formula we see that since both F, and G, are real, then 


: (—1) ote ne { =) = ; } 
G.= Pe Va No (22) r- 
L ira +1 cm R.I ( > Nor. s(—12) (A16) 





- 
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We can obtain the real part (R.P.) of (—ix/2)N2i',(—iz) from either (19) or (23) if 
we make use of the formulas“ 
(—iz/2)"J,,{(—tz) =e (x/2)"I,n(x), (A17) 
(—ix/2)"Y,,(—ix) = (—2/x cos rm)(z/2)"K,,(z) — i e~**"(2/2)"I,,(2). 
Using (23), we evidently arrive at the asymptotic expansion of G, given by Breit and 
Hull’, and from (19) we obtain the convergent series of Jackson and Blatt®. We have 
explicitly checked the terms containing (1/n*) to the powers zero and one in our resulting 
expansions for L = 0 against those given by Breit and Bouricius’ and Jackson and 
Blatt. The agreement is complete except that the latter authors have omitted a minus 
sign in the term H,(r) of their equation (A3.14); the sign of this term is given correctly 
by Breit and Bouricius. 


M4] bid., pp. 372-373 
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BOOK REVIEWS 


Abwickelbare Flachen. By Hans Schmidbauer. Springer-Verlag, Berlin-Géttingen- 
Heidelberg, 1955. V + 66 pp. $1.57. 
The profusely illustrated booklet treats the elementary descriptive geometry of developable sur- 


faces with particular attention to problems of engineering design. 
W. PRAGER 


Relativity: the special theory. By J. L. Synge. North Holland Publishing Co., Amsterdam, 
and Interscience Publishers, Inc., New York, 1956, xiv + 450 pp. $10.50. 


The special theory of relativity is now over fifty years old, and it is here treated as an established 
vecial pleading for its acceptance. Professor Synge spends no time on the usual 


theory that needs no sg] 
preliminary discussions of ether, absolute motion, Einstein’s operational definition of simultaneity, and 
the like. He begins with the idea that the physical world can be pictured in terms of a four-dimensional 
metrical continuum, and on the basis of quite mild assumptions about the temporal ordering of events 
shows that the local metric must be of the Minkowskian type, the Newtonian type appearing as a limiting 
case. He then makes the specific assumption that the space-time is flat and proceeds to develop the special 
theory deductively, pausing at frequent intervals to make contact with the familiar, though not wholly 
appropriate concepts of Newtonian physics, but returning always to the Olympian viewpoint of Min- 
kowski. The first three chapters set the Minkowskian stage. There follow two chapters on the Lorentz 
transformation, the first discussing it mathematically in considerable detail, and the second applying it 
to derive such familiar results as the contraction of lengths, and also to construct a special relativity 
model of an expanding universe. Subsequent chapters take up the mechanics of a particle and collision 
problems, the mechanics of discrete systems, the mechanics of a continuum, the electromagnetic field 
in vacuo, and fields and charges. There are five appendices discussing various technical points, a list of 
references, and an excellent index. 

Of the wealth of interesting items in this book only a few can be mentioned here. Thus the author 
emphasizes the relativistic sameness of photons, and comes to the logical, if at first startling, conclusion 


that luminosity should be defined in terms of all photons and not just those of visible light in a given 


frame. Again, he stresses the difficulties attending the definition of a rigid body in relativity, and avoids 
having to face them by postulating an ingenious mechanism for the transmission of impulses between 
particles. He shows how center of mass, a by no means obvious concept in a theory that lacks absolute 
simultaneity, can be defined with the aid of angular momentum. He shows how conservation of angular 
momentum can ameliorate an awkward and surprising lack of uniqueness in certain collision problems. 
He gives detailed discussions of many collisions involving material particles and photons that are of 


special interest to nuclear phy sicists, contrasting the various relative, three-dimensional pictures with 


their absolute, four-dimensional counterparts, though here he may confuse the reader by his failure to 
between space-time diagrams of world lines and those of vectors. He makes the treatment 


distinguish 
examining it in four different frames of reference. 


of the Compton effect particularly illuminating b 
Out of the pure electromagnetic field he constructs a remarkable non-singular model of a material particle 
which, because it is free of singularities, must unfortunately be uncharged. He makes every occasion an 
four-dimensional geometrical elucidation. And he does many other things of which the 


opportunity for 
above is a wholly inadequate sample. 
Despite its 450 pages, the book confines itself to what may 


from the quantum aspects of the special theory of relativity. Though spinors are not wholly absent, 


thev enter only briefly in connection with the Lorentz transformation. 
While the above may indicate the scope and nature of the book, it does not convey the author’s 


be called the classical as distinguished 


characteristic blend of insight and gentle wit that gives the book its special quality. The explanations 


are patient and penetrating, and remarkably lucid, often being expressed in the scientific equivalent of 
words of one syllable. In the preface the author writes that his aim is “to make space-time a real workshop 


for physicists, and not a museum visited occasionally with a feeling of awe.” In this he has admirably 


succeeded. 
BaNesH HOFFMANN 


(Continued on p. 64) 
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THE EFFECT OF TRANSVERSE SHEAR DEFORMATION ON THE 
BENDING OF ELASTIC SHELLS OF REVOLUTION* 


BY 
P. M. NAGHDI 
University of Michigan 


1. Introduction. ‘The classical theory of thin elastic shells of revolution with small 
axisymmetric displacements, due to H. Reissner [1] for spherical shells of uniform thick- 
ness, and Meissner [2, 3] for the general shells of revolution, was recently reconsidered 
by E. Reissner [4], where reference to the historical development of the subject may be 
found. Although the formulation of the linear theory and the resulting differential 
equations contained in [4] differ only slightly from those of H. Reissner and Meissner, 
they offer certain advantages not revealed in earlier formulations. 

As has been recently pointed out, the improvement of the linear theory of thin 
shells, by inclusion of the effects of both transverse shear deformation and normal stress, 
requires the formulation of suitable stress strain relations and appropriate boundary 
conditions which, for shells of uniform thickness, have been very recently carried out 
by E. Reissner [5] and the present author [6]**. The latter also contains explicit stress 
strain relations when only the effect of transverse shear deformation is fully accounted 
for but that of normal stress is neglected. 

The present paper is concerned with the small axisymmetric deformation of elastic 
shells of revolution, where only the effect of transverse shear deformation is retained. 
The basic equations which include the appropriate expression for the transverse (shear) 
stress resultant due to the variation in thickness, are reduced to two simultaneous 
second-order differential equations in two suitable dependent variables. These equations 
are then combined into a single complex differential equation which is valid for shells 
of uniform thickness, as well as for a large class of variable thickness. Finally, by an 
extension of the method of asymptotic integration due to Langer [8], the general solution 
of the complex differential equation is discussed. 

2. The basic equations in surface-of-revolution coordinates. ‘The parametric equa- 
tions of the middle surface of the shell may be written as 


r=r(), 2z =a), (2.1) 


where £, together with the polar angle @ in the z, y-plane constitute the coordinates of 
the middle surface. Denoting by ¢ the inclination of the tangent to the meridian of the 
shell, then 


dz 
tang = me (2.2) 
and 
r’ = a cos¢@, z’ =asin®g, (2.3) 


*Received January 20, 1956. The results presented in this paper were obtained in the course of 
research sponsored by the Office of Naval Research under Contract Nonr-1224(01), Project NR-064-408, 
with the University of Michigan. 

**The stress strain relations derived in both [5] and [6] were obtained by application of E. Reissner’s 
variational theorem [7]. References to earlier work on the subject appear in [5] and [6]. 
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where 


a = [(r’)? + (2’)?]'” (2.4) 


and prime denotes differentiation with respect to &. 
The square of the linear element for the triply orthogonal curvilinear coordinate 
system £, 0, and ¢ (measured along the outward normal to the middle surface) is 


ds’ = a'(1 - i de? + (i + i) de +d, - (2.5) 
E Le 
where 
R=-3, R=-—— (2.6) 
sin @ 


are the principal radii of the curvature of the middle surface. 
For axisymmetric deformation of shells of revolution, the displacements in the 
tangential and normal directions may be taken in the form: 


U, = ulé) + cB), W. = wif6), (2:7) 


where u and w denote the components of displacements at the middle surface, and 6 
is’the change of the slope of the normal to the middle surface of the shell. As pointed 
out in both [5] and [6], the approximation (2.7) for the displacements is consistent with 
the neglect of the effect of transverse normal stress in the stress strain relations for a 


thin shell. It is convenient to express the displacements u and w in terms of the corre- 


sponding components u, and w, (on the middle surface ¢ = 0) along the radial and 
axial directions, respectively. Thus, 
Uu, = UCOSd — Wsin ¢, WwW, = using + w cosd (2.8a) 
or alternatively, 
u =u, cosd + w, sin ¢, w= WwW, COSd — u, Sin d. (2.8b) 


The components of strain, as given in [5, 6], with the aid of (2.8) may be written in 


the form 








on et ew on 
me s °= 
gE r ’ r ] 
(2.9a) 
cos } sin @ 
ve = wi —— — u' — + B 
a 
B’ r’B 
= e ’ Ke = ee (2.9b) 


and the relevant compatibility equation, obtained from (2.9a) by elimination of u, is 
(r’es) = (res)’ + 2'w, (2.10a) 


where 

w= Ye — B. (2.10b) 
It may be noted that, upon the neglect of terms involving y;, (which represent the 
effect of transverse shear deformation), (2.9) and (2.10) reduce to the corresponding 
expressions of the classical] theory. 
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The stress differential equations of equilibrium, given in [4], may be written as 
(rPy)’ = —Taqy , 
(rPy)’ — aNe + ragn = 0, (2.11) 
(rM.)’ — a cos¢M, — ra(Py cos¢ — Pysin ¢) = 0, 


where P, and Py are “horizontal’’ and “vertical” stress resultants, respectively; ¢z 
and gy are the components of the applied load in the horizontal and vertical directions; 
M, and M, denote the stress couples; and the stress resultants N; and Ng , as well as the 
transverse stress resultant V (due to 7;;), are related to Py and Py by 


aN, =1r'Py + at ae 
aV = —2’Py,+1r’'P,, (2.12) 
aN, = (rP x)’ aa Tadu . 


We recall that the stress strain relations (as well as the stress differential equations 
of equilibrium) which are employed in the classical theory of shells (and plates) of 
variable thickness, are those for shells of uniform thickness, although the resulting 
differential equations take into account the effect of thickness variation. As the present 
analysis is concerned mainly with the effect of transverse shear deformation, it becomes 
necessary to modify slightly the available stress strain relations of uniform shells [5, 6] 
by incorporating the effect of thickness variation into the expression for the transverse 
shear stress 7;, . To this end, we replace in Eqs. (2.9) of Ref. [6] the expression corres- 


ponding to r;; by 


(0+ fee Sof ~ GH] ~~) 
+E). em 


which, together with o; , satisfies the stress boundary conditions on ¢ = + h/2, where 
the direction cosines of the outward normal to the middle surface are {—h’/2a, 0, 1}/ 
(1 + (h’/2a)*]*. While expression (2.13) will not change the differential equations of 
equilibrium (2.11), it does contribute to the stress strain relations, as may be seen from 
the variational equation employed in [5, 6]. However, in view of the neglect of the trans- 
verse normal stress, together with neglect of second-order corrections in h/R in comparison 
with first-order corrections, the resulting stress strain relations for a uniform shell (in 
the form of (3.6) of Ref. [6]), except for a modification in the transverse shear stress 
strain relation, remain unaltered. 
We close this section by recording the stress strain relations in a manner suitable 
for subsequent analysis. Thus, 
0 1 


1 , 
= (Ne — Ns) + wn} (6 + Xs), 


— 1 (r’ 
=F (No — »N,) — ak (c B + vB )|. (2.14a) 











4 P. M. NAGHDI [Vol. XV, No. 1 
; 6 iyp_ uh) M, 
o 5Gh ny a! 
aM, o| (5 +o 3) _ an(" — “)|. 
\ r P / 
: (2.14b) 
id fu, 
aM, = D ( 8 + vB’) + ad{ 
NS r 
where 
E Eh’ 
C Th G= —, D = ——z 
Eh, 2(1 + »v) 121 — Vv)’ 


E and v are Young’s modulus and Poisson’s ratio, respectively, and 
’ . 


he re ) ; be - ‘ (2.15) 
12(1 — v)’ : R, RJ a@\r %)}: — 
3. Differential equations of shells of revolution. With 8 and (rP,) as basic dependent 
variables, proper elimination between Eqs. (2.12), (2.14), (2.11), and (2.10) leads to the 


following two second-order differential equations: 
F » (rD/a)’ Ke’ ) 
” — bk? —_ 9)/ , 
B+ k wn) En(2n +7 js 
( : (ey (r’D | ( k’ \() 
— — kr — yp | + wkd 20’ . _ 
\G kx*) : ” “= D/e) vkX\ 20’ + k ; B 


, 


Z , r f —_ 
+ — (rP,,) ( ) — 1P.)! — »(rP,)"” 
(rD/e) (rPy) + cll [(r’P») v(rPy)’"] 


7 E o op i A") Jer 
[O20 4. Blo 4 (£) lerar} 
+ [22 4 GBled _ (ler, 
- [ Bo8r + emey — (F) foam} 
sear + ED oriy — [(C) + YD ora 
= tna — cB dor + [2ft) + (2) +2) fo 
[OC +00+00+O8} 6.2 


- E ae (2'/aC) [ern + y= (rPy)’ 


rT (r/aC) 


k 


(3.1) 


(r/aC)’ “| —_— 
| ee +9! (ragu) — (radu)! — G7qGy VE » 
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where k and \ are given by (2.15), and (rPy) = —Jfr aqy dt. With 6 and (rP,) known, 
the radial and axial displacements may be determined from the second of (2.14a) and 


the equation 


.=/ (2c? + r'w) dé (3.3) 


which, with the aid of (2.3), is deduced from (2.9a). 

When the effect of transverse shear deformation is neglected (by setting \ = ye = 0), 
(3.1) and (3.2) reduce to the corresponding equations of the classical theory given by 
E. Reissner [4], and which were first derived in a slightly different form by H. Reissner 
[1] and Meissner [2, 3]. 

Introduction of the function ¥ defined by 


y = ae (Py), m= [121 —»*)]'”, (3.4) 


into (3.1) and (3.2) and some rearrangement of terms results in* 


™ r/a) ey  2{ A” h’ 7 
pt | a — a? (er +3%) — ax +)|s 
nce. vn) B~ pact) 
U a (“) (a2 +8 8 Rae \h /* 


1 tr (r/ /x)’ oe hae h’ 
+ vk ny +{™ + (a) y= +6 ly 





(3.5) 
1} (r’\’ ‘i (r/a)’ ) 
(Ys e (+ oe 
_(*¥ oo M (gM It & c/a)") ly} 
() » h al h (3 h vr + d + (r/a) v 
= F, ] 
me (r, a)’ ’ h’ = oe 
ie E a) i h alle ne 
7 [r! ; (r’/a)’ 1211 + - ry a se (2 ‘a)’ zr) 
|(£) eae i 5 r 2 2F (r/a) 2r 
— (ear(E - 2 Wily +e +o a 
— vk’ fe” + E + ti +2 PS at = (i a - a — ayy |e (3.6) 
' r/a) 





“iC 048) ey 8 
+ 1% (“la = a’ yr Js} 


= F,, 








*In Eqs. (3.5) and (3.6), as well as in the subsequent analysis, the effect of thickness variation due 
to shear stress strain relation is placed in double-squared brackets, i.e., [[ _ J]. 
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where hy is the value of the thickness h at some reference section (say = &); F, denotes 
the right-hand side of (3.1); and F, and T are given by 


2" 2 ) 2’/aC)’ ieee 
F, m EE (1 — 1201 + ” | 4 rs. ¥ oh ite k ‘ary |oP.) 
ov 2 


~ EW r (r/aC) 


2! . ” r/aC’)’ > 
+y - (rPy)’ — E + a rk 90)] |raands- aan)’ 


r (r/al) 


—— A/V) 1201 + 2 | , 
.* af h 5m i (3.7b) 


ata Et L) «a; 


m \R; R, 


Since h/R < 1, then 


hence, in what follows, terms of the order O(kX*) will be neglected in comparison with 
unity. With 
X,=B+hk'*rY, XX. = wy — vk, (3.8a) 


from which 
B = (X, — vk'yX,) + Ok’), 





(3.8b) 
y = (X, + vk'?rX,) + Ok’), 
and with the aid of the relations 
ey a 
(kr) h’ 
iin Madi (3.9) 
(k*/*n)”” aster 9 Nh’ , hi’ 
‘ae ee dh ? 
(3.5) and (3.6) may be written as 
y ho 1/ hoe ” 
L(X;) + A) + vk “rg: + an(u's he |x. 
(3.10) 
_ yay (2r’ ME \ys =P 
- (22 + ® h me Fy 5 


L(X;) = (y — K)Xo 4 
~ 2a + (otro + 9 — 1) (es) 
uty( 1 + (vk "Ng + G2) ([- TI] \e fz X, (3.11) 
a, (or Ah’ nT lye = FB 
™ E (22 ~— 35 ) * urn |x: = F,, 
where 


uw ya [OX yak | y 


(r/c) 
(“VY _ ©)! yin “j(™) ri 
| (£) ”“(/2) ila \ 9 Malas mh | ?» 
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1+» (“) _ fa)’ _ L(y’ = L(e ~) fey ™ oe: Ee + 74 
g y r ” “(r/a) vy \r v \r / (r/a) d A L(r/a) vr |’ 
= ry 1 (“)’ (r’/a)’ , 12(1 + ») (“) L(e ~) Stay" 
92 2 t vy \r +? (r/a) + 5 r ai v (r/a) 
we Ee : +4 
d A L(r/a) vr |’ 
oA WW gM lo)’ _ 4+ 37 () 
“es h + h E h 3X " “+i (r/a) v r/}y’ 
—— h’’ h’ h’ (r/ a)’ d’ _2 —» (“)] 
a= -3N El gh 43 (r/e) + $y vy \r/ i} 
of (r’/a)’ 4 6+ +9) 
"G/e) + 5 \r/ | 
| a aes (ley ")] 
ea (Ze r/ |’ 


2 
am 


u f(é) = Raho’ 




















a 
| 


(3.12) 





and it is to be noted that in the last of (3.12), u is a constant and f(€) is independent 


of thickness h(é). 
Multiplying (3.11) by 76 (where 4 is a function of &), adding the resulting equation 
to (3.10), and observing that 


X{ = (X, + 16X,)’ — 16X; — 18'X, , i = (-1)”, (3.13) 


will result in a fairly intricate complex differential equation. By taking 6 in the form 
; eines , : . ho =] 
i + ( vk’ ‘Mg2 + G2) — ((- 1)( ss ") | 
¥ rr my 
=—g%~ (es 7) 
h -1 
+ {{ + vk! (gi, + ay(x'f he) | (3.14) 
1/2 r’ 2 ho ins 
| 1+ [eh Xg2 + G2) — [[- TI] eS 


-[3o - oer8) TP” 


i’ = 6” = 0, (3.15) 


the complex differential equation mentioned may be reduced to 


with the restriction that 
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1/2, \/ Qr! , 
L(X, + 18X;) + il ok" (oe 2r 19 r) ” (rn fox. + i8X,)’ 
(k’“*X) h 


vr 
=e Pete 3 rr! «, fe dia f a 

_ iau'p("2)| + (a “AM(g2 + G2) — ([ ry) a4 he) ee + 16X,) 

(3.16) 
= F, + i6F 
+ | (22 tans =h~ = a (i+ 38%) 4. arn |x: 
vr h 
The coefficient of (X, + 76X,)’ in the above equation (say B) is 

ae ee (e ay 6 fr... k) = ie 
B= 7a) 7 3 i + il r (ky) oo ee [(7]] (3.17) 


and since 


6 = exp {2 / Be dgp = = (z) 1? exp in | 5(vk'’*d) 
: (3.18) 
- itv f [won(S — F) + tg mah 
then, by means of the transformation 
Y = O(X, + 16X,) (3.19) 
we finally obtain 
Y" + [Pye + AW))Y — OX; = O(F, + 76F,) (3.20) 


which, in view of the presence of the function X; on the left-hand side, may be called 


the “quasi-normal”’ form of (3.13). In (3.20), the functions &, ¥*, and A are given by 


, 1/2. \ 2r’ 2 r’ ¢2\ h’ | ¢2 | 2ry7 ‘ 
@(£) = (vk ’*d) nit (1 + y+) Ul — eee (1+ 36) | + [[67]] (3.21) 
, l 
and 
aaa (_ -_ : Poe oa’ 
W'*(é) = 13 + (vt A(G2 + G) — (I rs! | 
“7 (3.22) 


4 1. ( 2, h yr 4 i antes 
3 inf 7) h }! 


2 r’\° (r a)’ ‘ (r/a) 6 + ») (= : 
A(é) = (" )- 2 (r/a) - + 1] & a | 5 =) 
hae is 3 (r/a)’ h’ rh’ 
= 2 -) = 2 (r/ a) h + 3» rh (3.23) 
‘"\ (r/a)’ 1d" , Xd (3 ys fey") 
+ tévk’ ‘2 (E)' + (: ) Gi oa Tk er "S60 


(: 

re 
ne 5h” ie 5 (r/a)’ 9 d’ a 15 a} 
2h + ~2(rja) 2 2h 
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4. Approximation of 6. Since 6, as given by (3.14) is fairly complicated, we impose 
a further restriction on the order of magnitude of (u’)™", i.e., 


(u*)"* = O(k'*D). (4.1) 
Thus, by (4.1), vk*\u"? = O(kA*) which, as in the previous section, may be neglected 
in comparison with unity. With this approximation, (3.14) simplifies considerably and 


reads as follows: 


b= po alert) +f Lo aber TP. ae 
= oy — Nes; gy — Nes; - 
It should be noted that the restriction (4.1) is physically plausible, as it holds true for 
many cases of shells of revolution, e.g., ellipsoidal, paraboloidal and toroidal. 


We now return to (3.14) and observe that condition (3.15) is fulfilled if 6 is a constant, 


and this may be achieved by proper choice of y and x. 
For shells of variable thickness, as seen from (4.2), condition (3.15) is satisfied 


provided (y — x) (u’fho/h) ~* is constant. Hence, by (3.12), 


(ee | 7 LE] - mk & - S09 (2), (4.3) 


L 
which resembles the corresponding equation of the classical theory, first given by 
Meissner [3] and derived in a different manner recently by Naghdi and DeSilva [9]; in 
(4.3), K is a constant and the effect of the transverse shear deformation is represented 


by the second term on the right-hand side. 
For shells of uniform thickness, « vanishes identically and y in general will not be 


constant. However, for numerous shell configurations, 


s=1—hin(uty) (4.4) 
2 h 
may be approximated to a constant, in which case the coefficient function W’ reads 
ww? = [1 + wk' rl f (4.5) 
and (4.3) may be reduced to 
h = (me r,)| - ee Bs A (4.6) 


It follows from (4.6) that y and 6 will in fact be constant if both R, and R, are constant. 
It is clear that this requirement is more restrictive than the corresponding result of 
the classical theory [9] where h = (mK/v)R; . 

5. Formal solution of Eq. (3.20) by asymptotic integration. We conclude the present 
paper by discussing formally the solution of the homogeneous differential equation 
associated with (3.20), namely 

Y" + [fu + AJY = 06X;, (5.1) 
where uy is a large parameter. 

It follows from the work of Langer [8] that corresponding to the homogeneous equa- 
tion associated with (5.1), that is, the equation 


Y” + (Pu? + AJY = 0, (5.2) 
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there exists a related differential equation whose solution is asymptotic with respect 
to » to the solution of (5.2) and that its domain of validity is dependent upon the character 
of the coefficient functions of Y. If — ranges over an interval 7; which includes a point 
t) at which (i) the function A may admit a pole of first or second order*, (ii) Y° may 
contain as a factor the quantity ( — &)*°, a being any real positive constant, and 
(iii) both A and W’ are analytic bounded elsewhere in J; , then the asymptotic solution 
mentioned will be valid in the entire interval J; including & . If, on the other hand, 
the behavior of A and WV’ at & does not meet the required specifications, then the solution 
will be valid in a sub-interval of J; . 
Thus, if we write A and WV’ in the form 

ap fe ga + A,(), (5.3a) 

VW’ = & — &) Mi, (5.3b) 


A = 


where A, and B, are constants, A, is analytic and bounded with respect to uw in J; , and 

9 « . . . . . . . . . . 
Wi is a non-vanishing single-valued analytic function in J; including & , then according 
to Langer [8] the functions 


\yi = (¢ — ns 8e/9— 1g 1/81 2 | J,(¢)} (5.4) 
PA ly) 
are the solutions of the related differential equation 
” 3 272 A, B, ’ — 
y+ lee t+ €& — &)" + €— &) + 2) ly = 0, (5.5) 


where & is analytic and bounded with respect to » in J; , J, and Y, are Bessel functions 
of the first and second kinds, and 


p=7, b= (1-44), 
. (5.6) 
e= iy | VW) do. 
Writing A, as 
A, = a) + A), (5.7) 
where A is also analytic and bounded with respect to u in J; , then (5.2) may be written as 
y” + | eev'@ + orn, + 7m b) + ag |v = —A()Y. (5.8) 


Since the left-hand side of (5.8) is identical with the left-hand side of the related 
differential equation (5.5), then, by the method of variation of parameters, there is 


obtained 


Yu, =Y¥it D*(E, 0)[A(n) Vu;(n)] dn; j= 1,2, (5.9a) 
# ko 


*In cases of ellipsoidal and paraboloidal shells of revolution of uniform thickness, A contains a pole 


of second order at — = 0. 
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where 
D*(E, n) = Tr'(ys » yo)[ysEy2(n) — yolE)yr(n)] (5.9b) 
and I’, denotes the Wronskian of y, and y, . 
The integral equation (5.9) by the familiar process of successive iteration leads 


formally to the relation 
Yu, = yl) + Dye, (5.10) 
n=1 


where 


E ‘ 
y@ = [| DG, DAG)y!(n) an, 
in (5.11) 
yi; &) = y.(é). 
The proof for the uniform convergence of >>%., y;” in J; is given by Langer [8] and 
he has further shown that y; is dominant in (5.10) and thus, y; is asymptotic with respect 


to uw to the solution Y,, of (5.2). 
tecalling that 6 as given by (4.2) is a complex constant, i.e., 





5 = 6, + 16, (5.12) 
then, by (3.19) and (3.18), 
41, [¥'_ ey 
xi = 5, gm E o r|, (5.13) 


where gm denotes “imaginary part of’. Treating (5.1) as a non-homogeneous differential 
equation whose homogeneous solution is given by (5.10), then again by variation of 


parameters, there results 


Y, = Yu, - i Dre, nt O(n) P(n) dm | Fi (n) — 2 (n) vin } dn; (5.14a) 
j= 1,2, 
where 
Dt(—, 0) = Tr'(Vu, » Yn.) (Vu.@ Vu.) — Yu, ¥u.(n)] (5.14b) 


and I, is the Wronskian of Yy,, and Yy, . 
As in (5.10), by successive iteration, the integral equation (5.14) can be written 


in the form 


¥; = Yu, — D &@h”, (5.18) 
where 
(n+1) . r (n) 
KOI"? = [DG NOCH) LXE)]i” dr, (5.16) 


[X2©];” = [X2@)]; , 


and X} is given by (5.13). 
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or 
to 


The above solution is formally valid if }>*., [x]{” is uniformly convergent in J, . 
On account of the interdependence on \ of the various functions involved in (5.16), the 
proof for the uniform convergence of the series in (5.15) for general shells of revolution 
appears to be difficult and requires further investigation. Nevertheless, from the com- 
parison of (5.11) and (5.16), it appears reasonable to except that >> [x];” has in general 
the same behavior as >> y;”, so that 

\ 


; l 
= j + Ol- (O. 
Y; IH; r (2) 


and hence, y; is asymptotic with respect to u to Y; . 

Finally, it may be noted that whenever the right-hand side of (5.1) vanishes iden- 
tically, as in the case of spherical shells where \ = 0, the differential equation is of the 
same form as the corresponding equation of the classical theory (Ref. [9]), although the 


qr 
~] 
—_ 


coefficient functions are not the same. 
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THE ACCURACY OF DIFFERENCE APPROXIMATIONS TO PLANE DIRICHLET 
PROBLEMS WITH PIECEWISE ANALYTIC BOUNDARY VALUES’ 


BY 
WOLFGANG R. WASOW 
University of California at Los Angeles 


1. Introduction. The standard method for appraising the difference between the 
exact solution u of a problem involving a partial differential equation and the solution 
U of an approximating finite difference problem is based on the expansion of u by Taylor’s 
formula up to terms of third or higher degree. The partial derivatives of u that enter 
the argument in this manner are unbounded near the boundary C of the region R where 
the problem is to be solved, unless C and the prescribed boundary values are very smooth. 

In computational practice the boundary and the prescribed boundary values are 
almost always piecewise analytic. At the corners of C and at the points of C, where the 
boundary values have a jump in the first or second derivative the higher derivatives 
of u become usually unbounded. Therefore there prevails the unsatisfactory situation 
that most known appraisals of the truncation error U — u in the numerical solution of 
boundary value problems are based on assumptions which are hardly ever satisfied in compu- 
tational practice. 

The present note is meant as a first step to overcome this difficulty. It will be shown 
that for the simplest finite difference approximation to Dirichlet’s problem for Laplace’s 
equation the order of magnitude of the truncation error is not affected by jumps in the 
first derivative of the boundary function. Of course, this is true only outside the immediate 
neighborhood of those discontinuities. 

It is frequently contended that, as the data of a mathematical problem of physical 
origin are by their very nature only approximate, the discontinuities in the derivatives 
of boundary values can be safely ignored in problems of this nature. However, there 
are many problems where these discontinuities are a very accurate model of physical 
reality, while their elimination by a smooth connecting are would change either u or U 
in a manner that cannot be guaranteed to be small without a further investigation like 
the one given in this paper. 

It is hoped that the method of this paper can be extended to more refined finite 
difference approximations and to other differential problems as well as to problems 
where the boundary C has corners. 

2. Green’s function for the difference equation. Let C be a simple closed analytic 
curve on which is defined a continuous function f(s) that is piecewise analytic. By this 
statement we mean that C can be described by two analytic functions x = 2x(s), y = y(s) 
of period /, regular for real s, and that there exists a finite number of values of s, say 
0 <8, <s,--- <s, <1, such that f(s) is continuous, periodic with period J, and regular 
analytic in every one of the intervals s, < s < 8s,,,;,» = 1,---,n,ands, <s<s, +l. 
Furthermore, we require that (dx/ds)*? + (dy/ds)’ ¥ 0. 

Denote the interior of C by R. Then there exists a unique function u(x, y) for which 


Au = 0 in R, u=f on C. (1) 
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The function u(z, y) is regular analytic on C, except possibly at the points S, (v = 1, 
-++ , n) of C that correspond to the parameter values s, . 


’ 


We denote by A, the finite difference operator defined by 
A,U = h[U(a + h, y) + U(x — h, y) + U(r, y + h) + Ula, y — h) — 4U (a, y)). 


Let R, be the set consisting of all net points P in R such that the four nearest neighbors 
of P in the grid lie in R + C, and denote by C;, the set of net points in R + C which 
do not belong to R, . On C, we prescribe a function f, = f,(P) whose value at a point 
P of C,, is equal to f(s) at some point P’ of C for which PP’ < h. 

The problem (1) is to be approximated by 


A,U =0 in R, , U =f, on C,. (2) 


It is known that the truncation error v = U — u is O(h) if 0°u/dz* and d°u/dy’* are 
uniformly bounded in R. This was first proved by Gerschgorin in [1] by means of the 
maximum principle for the operator A, , which states that a function U for which A,U = 0 
in a set of gridpoints cannot have an extreme value in this set unless it is a constant. 
Since the third derivatives of u diverge near the points S, of C, a stronger tool than the 
maximum principle is now needed. Our study of the truncation error is based on the 
asymptotic properties, for h — 0, of Green’s function for the operator A, in R, . 
For the definition of Green’s function we consider the problem 


A.V = ¢(z, y) in R, , V 0 on C. (3) 


Its solution is clearly a linear combination of the values of g(x, y) at the points of FP, , so 
that we may write 


V(P) =? > GAP, De(Q). (4) 


QeRr 

Here we have written ¢(Q) for ¢(é, n), etc. The factor h® has been extracted in antici- 
pation of a comparison with Green’s function for the differential problem. G,(P, Q) will 
be called Green’s function for the operator A, in R, . If (8) is interpreted as a system of 
N linear algebraic equations for the values of V(P) at the N points of R, , the N° values 
of h’G,(P, Q) form a matrix which is the inverse of the coefficient matrix in the system 
for V(P). Since the latter matrix is symmetric, so is G,(P, Q): 

G(P, Q) = GQ, P). (5) 

By applying (3) and (4) to the particular function ¢(Q) = 6(Q, Q’), where 


(0, Q#Q’ 


6(Q, Q’) = (6) 
ty, 0. = @’ 
it is seen that G,(P, Q) is the solution of the problem 
A,.eG,(P, Q) = h-“i(P, Q), ren. 
(7) 


G,(P, Q) = 0 , Pe. 


The subscript P in the symbol A,.p means that the operator is to be applied with respect 
to the variable P. The function G,(P, Q) is non-negative in R, . For, G,(P, Q) cannot 
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be a constant, because of (7). If it is negative anywhere in R, , there must be a point 
P = P,in R, where G,(P, Q), as function of P for fixed Q, has a negative minimum, while 
G,(P; ,Q) > G,(Po , Q) for at least one of the four nearest neighbors P, , (j = 1, --- , 4), 
of P, in the grid. However, the difference equation in (7) implies that 


Gi(Po , Q) > } dX GP; ,Q, 


and therefore G,(P, , Q) must also exceed at least one of the four values G,(P; , Q), which 
is a contradiction to the minimum property of P» . 

The asymptotic study of G,(P, Q) is based on certain results contained in [2]. It will 
be convenient to write Sin and Cos instead of the usual symbols sinh, cosh for the 
hyperbolic function in order to avoid confusion with the trigonometric functions of the 
mesh length h. We shall be concerned with the function 

2 f" 1 — cos(rAje'”” 


( == : dy 
x(a, 7) Jo sinh pu . 





where yu is the function of \ defined, for 0 < X < z, by the equation 
cos \ + cosh wp = 2 


and the condition lim u/A = 1. It was shown by McCrea and Whipple [2] that x(¢, r) 
has’ the following properties. 


l. x(o, 7) = x(7, 0), 

2. x(0,0) = 0, 

3. xe+1,)4+x6¢-—1,)+x,7+ 1) + x(e,7 — I) — 4x(¢, 7) 
(0, r+o #0 


_——— 
4. x(o,7) = log (0? + 7) +¢4+ 0(c°), uniformly in r, 
T 
as ¢o — . Here, c = (1/m) (log 8 + 2y), and y is Euler’s constant. 
Actually, McCrea and Whipple give only 0(c~") as the order of the remainder terms, 
but their own calculations show that the stronger result is true. For reasons of symmetry 


we replace the error term by 0 [1/(c” + r*)], which is permissible in view of property 1. 
It follows, then, that the function 


1 1 
yi(z, y) = 1,(2,%) + 2 togh—e 


has the properties 


1’. vi(z, y) = vy; x), 


2’. (0, 0) 


I 


] 
on log h- Cc; 
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0, x+y <0 
3’. Axy.(z, y) = ) 


-2 


as z=y=0, 


, , 1 2 2\1/2 ( h? ) 
4’, (zr, y) = = log (2? + y”)'? + 0i-——). 
w(t, ¥) = 5 log (@ + y’) aa 
The function —y,(z, y) is a discrete analog of a fundamental solution for Laplace’s 
equation. We now introduce the function 
H,(P, Q) = —y¥(x _ e ioe n); (8) 


where (xz, y), (&, 7) are the coordinates of the points P, Q, respectively. Then we have, 


by property 3’, 
A,.pH,(P, Q) = h-*a(P, Q). (9) 


In order to compare the asymptotic behavior of H,(P, Q) with that of G,(P, Q) we 
introduce the difference 


e(P, Q) = GAP, Q) — HCP, Q), (10) 
which, because of (7) and (9), satisfies the difference equation 
A,. pe, (P, Q) = 0, Pe R,, (11) 


and the boundary condition 
e(P,Q) = —H,(P, Q), Pet, . (12) 
We shall show that 
e(P, Q) = WP, Q) + O(A), forP eR, , Q = P, and Q not on C, (13) 


where ¥(P, Q) is, for P e R, the harmonic function of P with boundary values 
VP, Q) =~ log PQ, Pel. (14) 
aT 


To this end we expand A, py(P, Q) by Taylor’s formula and use the fact that ¥(P, Q) 
is, for Q not on C, a harmonic function of P in the closed domain R + C, thanks to the 
analyticity of the curve C (see [5], p. 187). Then we obtain 


A,.pV(P, Q) = O(h’), PeR,, Q not on C. (15) 
Also, by property 4’ and formulas (8) and (14), 
VP, Q) = —H,(P, Q) + O(h), PeC,, P#Q. (16) 


The last two formulas are valid uniformly in Q, if Q is bounded away from C by a positive 
distance independent of h and if PQ > ah’”’, where a is a positive constant independent 
of h. If we subtract (15) and (16) from (11) and (12) we find 


A, ple(P, Q) — WP, Q)] = O(n’), PeR,, Q ¥~ P, (17) 
e(P, Q) — WP, Q) = Oh), rec,, P#Q. (18) 


It is well known, and easy to prove by means of the maximum principle for the operator 
A, , (see [1]) that the Eqs. (17), (18) imply indeed (13). 
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Finally, we replace e,(P, Q) by its definition (10), and use once more property 4’. 
Then we see from (13) that 


‘Gu « >. log PQ + ViP,Q) +01), PeR,, QP, QnotonC. 


The first two terms in the right member together constitute Green’s function in R for 
Laplace’s operator A. Hence, 
G(P, Q) = G(P, Q) + O(h), forP eR, , QP, Q not on C. (19) 

We repeat that this relation is uniformly valid, if Q is bounded away from C by a positive 
distance independent of h, and if PQ > ah'”. 

3. The behavior of harmonic functions near corners of the boundary values. We 
consider first the special case that C is the unit circle and that u(z, y) = u,(z, y) is a 
function from the sequence defined recursively by 


Um(x, y) = Re F,,(2), z=a2+ wy, 
, (20) 
iaijnimti=28, P= -i f F,_( dt, teR+C. 
1 


If 


we take the branch of log (1 — z) defined by | a| < 2/2. Then F(z) is regular analytic 
in R + C except at z = 1, and 


F(z) = uo(x, y) + w(x, y) = —a + 7 log p. 
On the boundary, 


r— 0 


u(x, y) = f(0) = — 5 0 < 0 < 2z. 


Hence, fo(@) is continuous except for a jump of size at 06 = 0. For m > 0 the functions 

. . > ° i@ 
U,(x, y) and v,,(z, y) = Im F,,(z) are continuous in R + C. On C we have z = e” and 
therefore 


d’ , \ ’ f iv 
— F(z) = F-,(2)e oe z2eC, 


i.e. the boundary values f,,(0) of u,,(z, y) have the derivatives 


f,, (0) = f,,-,(0) cos v@ — v,,-,(z, y) sinv@ + --- , (x,y) eC, 


Jm 


where the dots indicate terms involving higher subscripts. This proves the continuity 
of f,.’ (8) for » < m. For v = m, the first righthand term has a jump of size (—1)"z at 
6 = 0. The second term equals log p sin v6, which has the limit 0 at 6 = 0. All the other 
terms are continuous, so that the jump of (—1)"x characterizes the behavior of f,,.”’ (6) 
at 6 = 0. 

The partial derivatives of u,,(a, y) of order y < m are continuous in R + C. This 
follows directly from the definition of u,,(x, y) in (20) and the Cauchy-Riemann equa- 
tions. The mth derivatives are multiples of u(x, y) or vo(2, y) and are therefore O(log p), 
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at worst. If vy > m, the derivatives are—except for numerical factors—the real or imagin- 
ary part of d’-”/dz’”” log (1 — z). Their order of magnitude is, thus, 0(p’™”). 

Now assume that the boundary function f(s) = f(@) of u(z, y) is regular analytic 
on the unit circle C except at 6 = 0, (i.e. z = 1), where it has a jump of size K in the 
pth derivative, while the derivatives of lower order are continuous. The function f(6) — 
(—1)’(K/x)f,(@) has then continuous derivatives on C up to order p inclusive. If it has 
a discontinuity in the (p + 1)st derivative, this can also be eliminated by addition of a 
suitable multiple of f,.,(@), and so forth. Hence, there exists a linear combination 


f*(0) = f(0) + cof (8) + asfys:(0) +--+ + aefnse(O) 


with continuous derivatives of order p + q and, at worst, a jump in the (p + g + 1)st 
derivative. The harmonic function with boundary values f{*(@) on C possesses then 
continuous derivatives up to order p + q, inclusive in R + C (see [3], p. 243]. Therefore 


a" “uy ou, "ile 
ox" o| ax ox 7 Ye + -" ox" | 
Pau | O(log p), q=0 
= 0 5p? = (21) 
Ox a 
(0(p 3 q>0 
and, similarly, 
arte, | O(log p), q = 0, 
ay ant ~ aad 
= (0p oon a> ©. 


a derivative of f(@) at any other point 6 = @ on the unit 
circle C can be analyzed by means of auxiliary harmonic functions obtained from u,,(z, y) 
by a rotation through the angle @ . If more than one point with discontinuities in the 
derivatives of f(@) occur they can be handled simultaneously by adding a sum of ap- 


A similar discontinuity of 


propriate compensating harmonic functions. Finally, if C is not the unit circle R can 
be changed into the interior of the unit circle by a conformal mapping. The mapping 
function and its inverse transform harmonic functions in R into harmonic functions 
inside the unit circle, and conversely. Since C is analytic the mapping function is analytic 
in C + R and heene the transformed boundary function is harmonic at all points of the 


unit circle except those corresponding to the points S; ,j7 = 1, --- ,n on C. Furthermore, 


the orders of magnitude of the harmonic function in R near the points S, are the same 
as those of the image function in the unit circle. Hence, we have proved that the solution 
u(x, y) of problem (1) is harmonic in R + C except at the boundary points S,; , where f (8) 
is not analytic. At those points the order of magnitude of the derivatives for approach in 
R + C is determined by the formulas (21) and (22). 
4. The truncation error for Dirichlet’s problem. 
for the approximate solution of problem (1) by means of problem (2) is the solution of 


The truncation errorv = U — u 


the problem 
Aw = —A,u , in FR, , (23) 


v=f,— U, on (’,.. (24) 
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By Taylor’s formula we have 


] 3° 9° 3° 
A,yu(x, y) = E ;u(x + 6h, y) — Z u(x — Oh, y) + =a u(x, y + Oh) 
6 | 0x OX Oy 


a° | 
— ay u(x, y — Oh (25) 


with 0 < @ < 1. We know from the preceding section that 6°/dx* u(Q) and 0°/dy* u(Q) 
are 0(QS;*) near S; . For simplicity we consider only the case that there is no more 
than one singular point S; = S on C, and set QS = p. The extension to a finite number 





of such points is trivial. 

From now on we must subject the choice of grids to the important restriction that 
the distance of the points S; on C from the nearest grid line be at least bh, 0 < b < 1, where 
b is independent of h. Then the right member of (25) does not exceed K,hp~ in R, where 
K, is a constant depending on f(s) and on b, not on h, i.e. 


| Ayu | < K,hp™’. (26) 


In order to arrive at an appraisal of v we represent v as the sum v = v, + v, of the 


solutions of the problems 
Ai; = —A,u, in R, ’ (27) 


y, = 0 . on C, 


and 


Aw, = 0 ‘ in 2, 
(28) 


v, = f, — U, on C;, . 
Discussing v, first we represent this function in the form (4) and make use of (26), 
obtaining 


lv(P) | < hK, >> NGAP, Q)p”. (29) 


QeRa 


If we extend the definition of G,(P, Q) in an appropriate manner from the points of 
R, + C, to all points of R + C, the sum in the right member can be replaced by an 
integral. To this end we associate with every square of the grid that lies entirely in 
R + C the value of G,(P, Q) as function of Q, at its lower left vertex. At the remaining 
points of R we define G,(P, Q) as being zero. Then 


> 1G,(P, Q)p? = I/ GAP, Q)p? dQ + O(n). (30) 
QcRr - 

Since G,(P, Q) is symmetric in P and Q, we can appraise the integral above by means 
of the asymptotic formula (19), provided P is restricted to some proper closed subdomain 
R’ of R. Let 

R=R, +R.+R;, 
where R, is the circular region PQ < ah'” about P, the domain R, is closed and satisfies 
R'’ CR, +R, C R, and, finally, R; = R — R, — R,. 
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We show first that 


/ G.(P, Q)p * dQ = Oh log h). (31) 
R, 
To see this we observe that the function | e,(P, Q) | assumes its maximum on C,, because 
of (11). In view of (8), (12) and the definition of y,(z, y) it follows that | e,(P, Q)| remains 
uniformly bounded, for Q e R’, as h — 0. Since H,(P, Q) = O(log h), in consequence of 
property 2’, we see from (10) that 


G(P, Q) = Ui log h), (32) 


‘if Q ¢ R’ or if P e« R’ (the latter because of the symmetry of Green’s function). This 
proves (31). 


Furthermore 
[| G,(P, Q)p7? dQ = I| G(P, Q)p? dQ + I| O(h)p~? dQ, (33) 
Ra Ra Rs 


by (19). The first integral in the right member exists, because p ‘ is bounded in R, , and 
the last integral is O(h). Thus 


[/ GP, Q)p2. dQ = O(1). (34) 
Ra 


In the remaining integral {fr, G,(P, Q)p ° dQ we make use of the fact that G,(P, Q) 
is, by definition, zero whenever Q is a point of a gridsquare that does not lie entirely in 
R + C. The point S lies in such a square; in fact, it has at least the distance bh from the 
edges of this square by virtue of the hypothesis introduced in the paragraph after formula 
(25). Therefore, R; may be replaced in the integration by a subregion R¥% that has at 
least the distance bh from S. Using again (19), it follows that 


[[ e@, Qe a@ = ff ae, @e* a9 + ff OMe” aa. 
Rs 2° oe 


The first integral in the right member is bounded, because G(P, Q) vanishes on the 
boundary and is there continuously differentiable so that 


G(P, Q) < Kap, Pet’, Qek;s . 

The last integral in the right member does not exceed 

b | O(h) | hp7 dQ, 

Rs 

which is bounded. Therefore 

I/ G,(P, Q)p * dQ = O()). (35) 

Rs 
If (31), (34) and (35) are inserted into (30), the inequality (29) is seen to imply that 


v(P) = O(A), rin. (36) 
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We now turn to the appraisal of v, . In order to make our analysis of Green’s function 
available for this problem, we use Green’s formula for the operator A, , which can be 
written 

h? ¥> (VA,U — UAV) +h > (VIr,U — UT,V) = 0, (37) 
Ri Ch 
(see [4], p. 36). Here U and V may be any two functions in the grid, and the operator 
I',U is defined as follows: Let Q be a point of C, and denote by Q; , (j = 1, --- » < 3) 
the gridpoints in FR, at distance h from Q, then 


r,U(Q) = h"'[U(Q,) + --- + U@,) — vU(Q)). (38) 
We apply (37) with V = v,(Q), U = G,(P, Q) and find, using (7) and (28), that 
v.(P) = —h y T,.ofG,(P; Q)|[F.(Q) — u(Q)]. (39) 
QeCr 


Now, by definition, f,(Q) = f(Q’) where Q’ lies on C, and QQ’ < h. In analogy with the 
hypothesis introduced after formula (25), we have to restrict the choice of Q’ somewhat by 
requiring that the distance from S to the points of the segment QQ’ be at least bh uniformly 
for Q on C,, . By the theorem of the mean 


f.(Q) — u(Q) = f(Q’) — u(Q) = i u(Q’’”) - QQ’, 


where 0/dc indicates differentiation in the direction from Q’ to Q and Q” is a point 
between Q and Q’. We proved in Sec. 3 that the first derivative of u at the point Q” in 
any direction is O(log p’’), the letter p’’ designating the distance SQ”. With the help of 
our restriction on Q’ we shall show first that 


| f(Q) — u(Q) | < Ksh(| log p | + 1), (40) 


K, being independent of A and p. By assumption, p > bh and QQ” < QQ’ < h, whence 
QQ” < p/b. Therefore, p’’ < p + QQ” < p(1 + 1/b). Using our other assumption that 
p”’ > bh we show analogously that p < p’”’ + QQ” < p’’(1 + 1/b). The ensuing double 
inequality 


1 
log p — log (1 + t) < log p’’ < log p + log (1 r 1) 
) 
implies that 


l | 
log po” | < max { log p+ log (1 + 1) , | log = log (1 + 1) } 
| ! 


lA 


| log p| + log (1 +4), 
and, therefore, 


| 2 uQ’) OG | = (0 log 6”) < Ksh(| log » | + 
On C, the relation (19) can be used, if P e R’, and this implies that G,(P, Q) is O(h), 
for Qe C, . Using (38) and (40) we obtain then from (39) the inequality 
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v(P) | < Kh } h(| log p| + 1), (K, a constant). (41) 
Cra 


It is plausible that 
do (| log p | + 1) = O01). (42) 


The simple proof of this statement will be postponed to the end of this section. On the 
basis of (41) and (42) we have 


|v(P)| = Oh), PeR’. 


This relation, together with (36) completes the proof of our main result, which we now 
state as a formal theorem. 

Theorem: The truncation error v(P) = U(P) — u(P) corresponding to the approximate 
solution of problem (1) by means of the equations (2) is of order O(h), provided the following 
conditions are satisfied: 

(a) the boundary C is a simple closed analytic curve; 

(b) the boundary function f(s) ts continuous and piecewise analytic; 

(c) the distance from the singularities of f(s) to the nearest point on a grid line is not 

less than bh, (0 < b < 1, independent of h). 

(d) If f,(Q) = f(Q’), Q e C, then the distance of the singularities of f(s) from the 

segment QQ’ is not less than bh. 

The truncation error v(P) has the order O(h) uniformly in every closed subdomain of R. 


Proof of formula (42). We show first that the total number M of the gridpoints in 
C, is O(h~'). The analytic curve C possesses only a finite number of points where either 
dz/ds or dy/ds vanishes. These points divide C into a finite number of arcs each of which 
does not intersect the same grid line twice. Hence, if C” is one such arc and if the lengths 

) £ 

of its projections on the axes are L, and L, , respectively, then the arc C° possesses at 
most (L, + L, + 2)h™* points of intersection with lines of the grid. It follows that the 
total number of intersections of C with lines of the grid is O(h~"). Now, every point of 
C,, is an end point of a mesh side of length A that has a point in common with C. Hence, 
M = O(h""), as claimed, say, 


M < Lh". 


For sufficiently small h no closed segment of length h joining two gridpoints will have 
more than one point in common with C. Then exactly one of the two end points of such 
a segment belongs to C, . We now measure the arc length s on C from S and order the 
points P, of C, in such a way that r, < r, if and only if P,, is an end point of a segment 
that meets C at a point with smaller value of s than any grid segment ending at P,, . 
Then SP, < 2'” rh and, therefore, 


M 
> A(| log p | + 1) < DY AC| log 2'rh | + 1) 
Ca r=1 


lA 


L+h 
I (| log 2'7¢| +1) dt< Ks,  (K,;a constant), 
0 


which proves (42). 
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BOOK REVIEWS 
(Continued from p. 40) 


Boundary layer effects in aerodynamics. Proceedings of a Symposium held at the National 
Physical Laboratory on 31st March and Ist April, 1955. Her Majesty’s Stationary 
Office, London, 1955. ix + 400 pp. $4.00. 

This book is a collection of nine papers presented at an international symposium on boundary 
layers held at the National Physical Laboratory in 1955. Included in the compilation is an introductory 
address by L. Howarth and a summary of the discussion which followed each presentation. The papers 
covered most of the topics currently under investigation in the field, and their character varied from 
those of an applied mathematical and fundamental physical nature, to those concerned more than with 
the practical engineering aspects of boundary layer studies. This fine collection would seem to make a 
basic contribution at each of these levels. An interesting feature of the compilation proves to be the dis- 
cussions which appear in summarized form. Unfortunately, some of the utility of the book is lost by the 
fact that nearly all the papers have recently appeared in the journals, or in other collected aerodynamic 
studies. 

The following is a list of the authors and the general topic of the paper: R. Timman (three-di- 

mensional boundary layers), M. B. Glauert and M. J. Lighthill (axisymmetric boundary layers), M. 

Gregory, J. T. Stuart and W.S. Walker (stability of three-dimensional boundary layers), G. B. Schubauer 

and P. 8S. Klebanoff (transition), D. Kiichemann (viscosity effects on swept wings), R. C. Pankhurst 

(boundary layer control), A. D. Young and 8. Kirby (profile drag at supersonic speeds), D. W. Holder 

and G. E. Gadd (shock-boundary layer interaction), H. H. Pearcey (transonic turbulent separation on 


airfoils). 
RonaLD F. PROBSTEIN 


Thermodynamics and statistical mechanics. By Arnold Sommerfeld. Academic Press, 

Inc., New York, 1956. xviii + 401 pp. $7.00. 

For anyone who is familiar with the Lectures on Theoretical Physics—Vols. I, II, III, IV, VI, by 
A. Sommerfeld—it is a sufficient review of the present book merely to say that it is Vol. V of the series. 
The remarkable understanding of theoretical physics by the author is in this volume as in the others 
made available to the reader by the clear and simple way in which the ideas are ordered and presented. 

As the name implies the book presents both the macroscopic thermodynamics and the classical 
and quantum statistical mechanics. About the first 2/5ths of this 400 page book is on thermodynamics. 
The last 3/5ths is on statistical mechanics. At the end of the first portion is a brief presentation (16 
pages) of “Irreversible Processes: Thermodynamics of near-equilibrium processes’’. 

Any book of 400 pages on a subject with an extensive literature must be selective in its material. 
The reviewer regards the selection of material for this volume as good. The first chapter on ‘“Thermo- 
dynamics—General Considerations” is followed by Chapter II on “‘The Application of Thermodynamics 
to Special Systems”. In the closing section of this chapter irreversible processes are treated. “The Ele- 
mentary Kinetic Theory of Gases” is treated in Chapter III. Next follows a long Chapter IV on ‘‘General 
Statistical Mechanics: Combinatorial Method”. Here the basic ideas and their applications to a number 
of special problems is treated in most lucid detail. The last, Chapter V, gives an “Outline of an Exact 
Kinetic Theory of Gases’’, where applications are carried through Conductivity and the Wiederman- 
Franz Law where the author’s own basic contributions are presented in proper perspective. 

As many of the readers of this review will know, the author died before he could finish this book. 
Several sections were completed by Bopp and Meixner who are to be congratulated on their ability to 
finish a work by a master book writer and great theoretical physicist without discontinuities showing. 
Although the reviewer has not attempted a comparison of the original German and translated English 
versions, the latter reads exceptionally well and gives every indication of being a faithful as well as 


fluent translation. 
H. Emmons 


(Continued on p. 82) 
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BIHARMONIC EIGENVALUE PROBLEM OF THE SEMI-INFINITE STRIP* 


BY 
GABRIEL HORVAY 
General Electric Research Laboratory, Schenectady, New York 


Abstract. A basic problem in the evaluation of residual stresses in simple elastic 
structures concerns the determination of the stress and deformation state produced 
by self-equilibrating, but otherwise arbitrary, normal and shear tractions acting on 
the edge x = 0 of a semi-infinite elastic strip (0 < x < », —1 < y < 1) which is free 
along the edges y = +1. This strip is known to experience, in accordance with St. 
Venant’s principle, inappreciable stresses at distances x 2 2 from the loaded edge, in 
spite of the very large stresses it may experience in the vicinity of the edge. An earlier 
paper [The end problem of rectangular strips, J. Appl. Mech. (1953)] based on the vari- 
ational principle, established approximate eigenfunctions (modes of response) and 
eigenvalues (laws of oscillation and decay) for the various possible self-equilibrating 
end tractions. In this paper we give a rigorous solution of the end problem. This solution 
is obtained in two steps. First we solve the two ‘mixed’ end problems: the parallel 
edges y = +1 of the strip are free, and along the vertical edge x = 0 (a) the shear 
displacement is given, the normal stress is zero, (b) the normal displacement is given, 
the shear stress is zero. These two problems are solved by extending the strip to the left, 
to —o, and finding the tractions that must be applied at y = +1 (x < 0) and at 
x = —,s0 that one have o, = 0,7 = 0, respectively, at x = 0, while the edge values 
of the displacements (more specifically, of dv/dy and u) are orthogonal polynomials 
in y (Horvay-Spiess polynomials and Legendre polynomials, respectively). The corre- 
sponding stress functions K,(x, y), J,(z, y) are found in the form of Fourier integrals 
plus polynomial terms. For x > 0 they may be rewritten as real parts of } ik Cre Px ; 
pm C,,, ®, , where ®, = z,’e*"* (cos 3y — y cot z sin z,y) or 2, e ““* (sin zy — y tan 
z, COS Z,y), and sin 2z, + 22, = 0. An alternate procedure for determining the coefficients 
Crear» Coax , based on a formula of R. C. T. Smith, which bypasses the extension of the 
strip to x = —~©, is also furnished. The second phase of the solution of the “pure” 
end problem—along the short edge (a) the shear stress is given, normal stress is zero, 
(b) the normal stress is given, shear stress is zero—consists in recombining the biharmonic 
eigenfunctions K, , J, within each class into functions H,(z, y), G,(z, y), so that the 
x = 0 values of H,, ., , G,.,, constitute two complete orthonormal sets of (transcendental) 
functions in y into which the given boundary stresses may be expanded. 

1. Introduction. We propose to solve the biharmonic eigenvalue problem of the 
semi-infinite strip. More specifically, we shall establish functions H,(z, y), G,(2, y) 
(even in y for even n, odd in y for odd n) such that 


(a) H, , G, are biharmonic functions 


V‘H, = V'G, = 0. (1) 


*Received January 27, 1956. 
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(b) H, , G, satisfy homogeneous boundary conditions of zero normal stress, zero shear 
1 
stress along the long edges (star denotes boundary value at y = +1)’: 


Hic2 = Ghz = 90, Adin = Gia, = 0. (2) 
(c) H, gives zero normal stress, G, gives zero shear stress along the short edge (the 
“degree” sign denotes boundary value at x = 0) 


H’ ,, = 0, Ge ., = 6. (3) 


(d) The edge values 
ty) = —Heev, SY) = Griv (4) 
constitute two complete orthonormal sets of functions into which prescribed self-equi- 
librating edge tractions 7°, o% , i.e., tractions satisfying the conditions 
+1 atl 


. +1 
[1° dy = [ w= / o°y dy = 0 (5) 
1 —1 =} 


may be expanded: 


r°(y) = D> (7°, teta(y), oy) = Do Co? , 8,)8a(y). (6) 


It follows that 


H(z, y) = > (r°, t,)H,(x, y) 
" (7a,b) 


G(x, y) = >> (2 , 8,)G,(x, y) 
are the stress functions of the two problems. We used the notation 
+1 
Go=f[ foowd, Islaan. (8) 
-1 


We arrive at the solutions H, , G, of the “pure” end problems by first solving the 
easier, “‘mixed”’ end problems pertaining to determination of stress functions K, and 
J,, which comply with conditions (a), (b), and with the modified conditions (c’) and 
(d’)’: 

(c’) Ke vy — 0, Ie zz - 0, (I, z = Ja); (9) 


si K°.. ~ y?, y = y’, ee 8” + 7y'*, ins (n of 2,3, 4, iad (10) 


ll 


ee (—1 + 3y’)/2, (—3y + 5y")/2, ee (11) 


The functions K°,, will be recognized as the (unnormalized) Horvay-Spiess poly- 
nomials Q, , Q; , Q, , --- (denoted formerly, except for normalization factors, by fy , 


1Clearly, because of the evenness or oddness of the functions H, , G, , the data at y = +1 also 
specify the data at y = —1. We use as distance unit the semiwidth of the strip, as stress unit the modulus 


of elasticity. 
*The advantage of regarding J, rather than J, = J,, z as the basic function will become apparent 


in the sequel. 
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fi, fe, +++ 3 see [1], Eq. (12)), with leading coefficient chosen as 1, which constitute a 
complete orthogonal set with respect to the boundary condition 

Q.(+1) = 0 (12) 


while the functions J2,, are the well-known Legendre polynomials P, , P; , Py, -:- 
forming a complete orthogonal set with respect to the boundary condition 


P(+1) = 1. (13) 

We disregard the first two Legendre polynomials, Pp) = 1, P,; = y, which represent 
rigid body motions, and do not give rise to stresses. In contrast, the singular functions 

Qo _ i, Q: => (14) 


which violate condition (12) (the condition of ¢,(0, +1) = o,(0, +1) = 0) are of con- 
siderable interest. It should be remembered that the corresponding stress functions, 
K, , K, may be written as linear combinations >>%3 c,K, of the complete set of functions 
K, , n > 2, and so their separate consideration is somewhat redundant. Nevertheless, 
a direct determination of Ky , K, is of great practical value; we shall therefore list their 
direct formulas along with those of K, , K,, --- . (The symbol Ky was previously de- 
noted by 2® in [4] and by 23, in [7], the symbol J, was previously denoted by 66 in [7].) 

It is clear that solution of the (a), (b), (c’), (d’) problem resolves the end problem 
of given shear displacement, zero normal stress, and given normal displacement, zero 
shear stress. For, if K is a stress function, then 


) - Oz ’ | ad VK vy + dv/dy (15) 
and, because of condition (c’), the edge values 
a) o/ 
| _ |“ og (16) 
Kew 0 
are properly specified. Similarly for a stress function J, we have 
Te=—f dy, Ty =e +4, (17) 
-1 


hence, because of condition (c’), the edge values 


| . | ” 
1 bot ue 
are again properly specified. Thus, the two mixed end problems have the stress functions 
, — (Kp2z » dv°/dy) ,- 
K(z,y) = >> (Ke. : L u) K,(2, y) 
2 | Kees | 


= (I> yy» U°) 
J(z,y) = 1,{z, y) = eww wt fo 
(x,y) = 1x, y > ie, (z, y) 





(19) 


as solutions’. 


‘It is obvious that the edge value problem—I and \/*J are specified along z = O— is also solved 
in terms of Eqs. (19). 
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Once we have determined the functions K,,(z, y), J,(z, y) appropriate to the mixed 
end problems, it is easy to calculate the shear stress boundary values —K°,, and the 


normal stress boundary values J°,, . We shall obtain these in integral form 


ae 2 ch dry sh ry , dy a 
7. = A. = i | aon od + @(A)2 ch 4 + P(A, | x? (20) 


of = J°.. = similar, 
where @(A, y) represents polynomial terms in \ and y, progressing to such power in \ 
as to compensate for the singularities introduced by the terms @/X", ®/X” (see Table IT). 
Knowing the integral representations of r9(y) it is easy to orthonormalize these expres- 
sions, writing 

to(y) = beerS(y), 

ts(y) = bast S(y) + Ba2rQly), (21) 


t.(y) best 3(y) + best S(y) + best S(y) 


and choosing the b,, so that (6,,, = Kronecker delta) 
‘if = &, (22) 
It follows that 
Hz, y) = b.K,.(z, y), 
H(z, y) = bu K,(x, y) + by2K.(x, y) (23) 


are the stress functions which resolve the edge value problem r° = given, o2 = 0, in 
the form of (7a), and 

G.(z, y) = Q22J.(z, y), 

G(x, y) = AgJd s(x, Y) + Ay2J2(z, y), (24) 


where the a,, are determined from the requirement that 


s.(y) = A220 3(Y), 


II 


8,(y) Ayo9(y) + Ay2a3(y), (25) 


(8, ° Sm) = Oni 


are the stress functions which resolve the edge value problem oc = given, r° = 0, in 


the form of (7b). 
2. The functions K,, , /, , H, , G, . Let z, denote the first quadrant roots z, , 2, --- , 


and 23 , 23, °°° Of 


sin 2z, + 22, = 0, sin 2z, — 2z, = 0 (26) 


1957] BIHARMONIC EIGENVALUE PROBLEM OF THE SEMI-INFINITE STRIP 69 


respectively. We recall that the Fadle-Papkovitch solutions [5], [6], [2], of V‘® = 0 
(z > 0), 
2, 4, 6, +++) 


(x, y) = 2, e “"(coszy — y cot a sin xy) (k 


(27) 


®,(z, y) =z, e "(sinzy — y tan z, cos zy) (k = 3, 5,7, ---) 


satisfy the homogeneous stress conditions (2) along the long edges of the strip, but 
produce both «2 and r° values, and both u° and v° displacements along the short edge. 
What is desired are, however, such combinations of ©, which give either o° stress or r° 
stress, making the other stress zero*. Such are the functions K, , J, and H,, G, . 

In the remainder of this section we give the results of the analysis..The analysis 
itself will be carried out in Secs. 3 and 4. 





In Table I below we list, up to n = 6, the expansion coefficients in terms of the 
Fadle-Papkovitch functions ®, , of the eigenfunctions 
K,(2, y) = ® DI Cut(z, y) (28) 
7 
TABLE I. Formulas for the coefficients Cyx in the expansions of K, and I,, , Eqs. (28), (29). For n = 
even, C!, Cr cos*z,/sin2z, , for n = odd, Cl, = Cix/cos 2 is shown. The subscript k of z% is omitted 
for the sake of simpler notation. 
K,, dv° dy Cee 
Ky, | 1 2 
K,\y | —2e7" 
re 1 — y 4z°°(2 + cos’ z) 


K,|y-y | —42°*(6 — cos’ z) 
‘ 


K,|1— 8y’ + Ty 82" *[42(3 + 2 cos’ z) — 2°(34 + 3 cos’ 2)] 
K; | y — 4y* + 3y’° | —82 *[30(9 — 2 cos’ z) — 2°(18 — cos” z)] 
K, | 1 — 19y? + 5ly' — 33y’ | 322°°[1485(4 + 3 cos’ z) 


— 1827(111 + 19 cos’ z) + 22°(13 + cos’ z)] 














a Th | ec. 
I, | (—1 + 3y*)/2 | —62"? 
I, | (—3y + 5y’)/2 3027° 
I, | (3 — 30y’ + 35y*)/8 30z~*[7(2 + cos’ z) — 32*] 
I, | (15y — 70y* + 63y°)/8 —21027*[3(6 — cos’ z) — 2°] 
2 | (—5 + 105y? — 315y* + 231y°)/16 | —210z2~°[99(3 + 2 cos” z) 
— 62°(15 + 2 cos’ z) + 22*] 





4A similar problem may be formulated also with respect to the end displacements u°, v°. See, how- 


ever, the footnote". 





70 GABRIEL HORVAY [Vol. XV, No. 1 


appropriate to dv°/dy = given by (10), c2 = 0, and of the eigenfunctions 


I(x, y) = ® 2 Cu®,(c, y) (29) 
k 
appropriate to u° = given by (11), r° = 0, where 
R = “real part of” (30) 
a for n = even, a for n = odd 


S9.f4.°** 


TABLE II. Formulas for the edge values of the stresses. 











K, ?. 
| { ao 
Kp —=| k,. ad 
T Jo 
g.i —§ [theo + kro] dX/ 
* ve 
| 
| gr ‘ 
i = | (2k,. — kre] dd/X? 
j T Jo 
_1:ser. — 2 
K; = [ [6keo + (6 + X°)kyo] dA/d 
= 
oS ial (4 + ane. + Env — 9) | ar 
T “& 
- 16 2 97 Oy 2 4) 7, 3 4 2 4 J, 6 
eo} woe ef (270 + 18N*)keo + (270 + 78° + A')kwo — FAL — 6y* + 5y') | dd/d 


K, 64 [ | ( 5940 + 1998” + 26r*)k,. — (4455 + 342” + 2a‘)k,, 
T 


+ 2. y(1 — y495 — aKe - 33y/)} | dx/d° 











Fe oS 
Js | [ jee Od/r* 

1, | 2 P Daeo + deol AN/X' 

Js | 80 [ (14 + 3d*)j,. — Tre] dd/r* 


| ® Jo 


= 


4? ai 
AO [08 + dijo + (18 + 4N°)je0] A/D’ 





- _ 420 [ [ (297 + 90d? + 2A*)j,. — (198 + 12d7)j,, + 33 (l- ay) | dr/r° 
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and summation extends over the first quadrant roots of (26). The two procedures, one 
due to the author, the other based on a method of R. C. T. Smith, by which the co- 
efficients C,,, may be determined, are described in Secs. 3 and 4, respectively. 
Anticipating the results of Sec. 3, and utilizing the functions k(A, y), j(A, y) listed in 
Appendix I, the edge stresses r° appropriate to K,, and the edge stresses o° appropriate 
to J, are found to have the integral representations listed in Table II. We give these 
expressions up to n = 6. Numerical integration, performed as in [4], then leads to the 
edge values of the stresses 73, 7 , --- , o4 displayed in Table III. The scalar products 
(r°r9), (¢2e°%), obtained by numerical integration over the tabulated values, are listed 


in Table IVa’. 


TABLE III. Edge values of the stresses. 








y= 0 .2 4 6 8 9 .95 io 
—— —— 

T) 0 —.10038 —.20638 —.32458 —.46371 —.54490 —.58929 — .63662 

TY . 18917 . 17263 .11981 .01867 —.16549 —.32227 —.43786 — .63662 

7) | 0 — .14339 —.27328 —.36970 —.38606 —.31681 —.22945 0 

73 . 15800 .13177 .05739 —.05029 —.15417 —.16970 —.14317 0 

7; 0 — .3644 —.5143 —.29678 . 27376 .54880 . 56353 0 

o> . 58820 . 53678 . 37249 .05747  —.51596 —1.0022 —1.3588 —1.9665 

o3 0 .47538 .81663 . 83568 .1329 —.86634 —1.7833 —3.7827 

o% — .87136 —.57794 .19518 1.0672 1.0413 —.25672 —1.8781 —6.3727 


TABLE IV. (a) Scalar products of 7°, 0%. (b) The expansion coefficients baz , Ant Of tn , Sn in terms of 
tT; , 0, , Eqs. (21), (25)*. 





nk (rorp) (oSo;) | Onk Onk 

22 . 1598 . 7912 2.501 —1.124 
33 .02726 1.590 6.056 — .7931 
44 | . 2848 2.722 | —1.939 — .6279 
42 | .05470 3836 | .6634 . 3044 


The orthogonalization processes (22), (25), finally lead, in accordance with (23), 
(24), to the expansion coefficients b,, , a,, of H, , G, in terms of K, , J, , as shown in 
Table IVb. The boundary stresses ¢, , t; , 44 produced by H, and the boundary stresses 
8S , 83 , §, produced by G, are plotted in Figs. 1 and 2. For comparison, the corresponding 
distributions based on orthonormalized self-equilibrating polynomials are also shown. 
(It will be recalled that the f,(y) polynomials are orthogonal, the fi(y), fi’(y) polynomials 
are not. They may be readily orthonormalized into functions 


7 = p Burfi ’ S, = > An i’ (31) 
2 2 


as shown in [2].) Note that the disagreement of the two sets of curves is not a measure 
of the inaccuracy of the variational approach, but reflects merely a rotation in function 
space from one set of orthogonal axes to another. 

*The last two digits of the entries of Tables IVa,b are uncertain. See footnote‘. 

‘Because this numerical integration was carried out over the sparsely determined values, at y = 0, 
.2, .4, .6, .8, .9, .95, 1.0, the entries of the present Tables IVa,b must be regarded as preliminary esti- 
mates. However, a more accurate determination would require a tremendous investment in time and 
personnel; these are not available to the author. 
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Fic. 1. The first three orthonormal boundary shear tractions. Full lines: 


t(y) = —H2.y = — >, by:K2.a ; 
i=2 


> Bu ift . 


dashed lines: T’.(y) 





[" 


$3 








S2 
Se 


S3 
So 
































“1.0 


Fic. 2. The first three orthonormal boundary normal tractions. Full lines: 


n 
8,(y) = Ge vy = } * Anil ? yy } 
i=2 


dashed lines: S,(y) = pe A,,f?’. 
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3. The method of analytic continuation. This method of approach was discovered 
in a somewhat accidental manner. In [4] we investigated the problem of the stress state 
in an infinite strip —~ <2 < +o, —1 < y < +1, occasioned by the temperature 


distribution 


a for «<0 (32) 
0 for x>0 
or, what amounts essentially to the same thing, by the edge tractions 
fe + , z <0 (33) 
0 z>0 


and found, incidental to the solution of the thermal stress problem, also the solution 
+ Ko(a, y) of the semi-infinite strip for the mixed edge conditions 


o° = 0, dv°/dy = 1/2. (34) 


Similarly, solution in [7] of the thermal stress problem of the infinite strip for 


0 z>0 
or, what amounts essentially to the same thing, for the traction distribution 
ra 
* = 0, ot "=? (36) 
0 zr>0 


provided also solution J,/6 of the semi-infinite strip problem for the mixed edge con- 
ditions 


re =0, wu? = (By — 1)/12. (37) 


It was then immediately obvious that all end problems of the semi-infinite strip 
should be reducible to problems of the (doubly) infinite strip. But instead of seeking 
determination, from a pair of integral equations, of complicated unknown tractions 
o* , r* applied to the left half, z < 0, of the horizontal edges of an infinite strip, in terms 
of the given values at x = 0 (dv°/dy = given, o2 = 0; or u° = given, r° = 0), one should 
be able to arrive at the tractions o*% , r* by inspection. This is the very program that was 
carried out in [4], [7] and [2], and led to the determination of the stress functions K, , 
J,, and K, , J; . The singularities of these four functions at x = — © are, however, 
not quite severe enough to illustrate fully the general approach. For this reason we 
determine below the stress function K, of the (doubly) infinite strip, appropriate to the 


boundary conditions 
P a > ° 2 S o- 
for 2 > 0: o, = 7 = 0 (38) 
for x=0: o2=0, adv°/dy = 1 — 8y’? + 7y’*. 

We shall find that K, has the form 


63 


5 Kuo — 42K, — 34K 2, — 6K2, , (39) 


K, = 
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where the functions K,, , --- , Ky, of x and y are defined, in terms of 
a= } = —dsh : 
A.A, = shv + dAchr, A.B, Ashi, (40) 
A.C, = shy, A,D, = —chi, 


(see the Appendix for A,) as follows:° 


94 p@ 
K,, = -— | | Ackay + Byshryy — 4 gin dz dy 
T Jo 2 , 
(41a) 
a i? a) a 
+4e-2/ (sin xx — dz 4% ref 
3 ad * 7 
K., = 2. [ | Coby + Dyshrxy — ly — | — dx 
T Jo 4 A 
(41b) 
> ig 
+(l1-y ae + =f (sin Ax — Az) al, 
P41 T Jo 
6 a@ 4 
K,, = J | | Ach\y + Byshrxy — ; a s (1 — y’) + [sine dy 
wT Jo 2 48 
_ 9 2 3.3 5 
+ 1642. +- | (sin Ax — Ar + = -- ve) al (41c) 
\6! T Jo 5!/7 AX 


2. 2 y) i se dr 
— =( (1 — y’) e+? += (sin Ax — Az) a), 





o> pe " 4 
x. «= | | Conny + Dyshrxy — * (1 -—y¥) + x (l—-y) ve. 


~* 9 a 
— 41 — is - 3 [ (sin xs =) a (41d) 


2 ier, Be dr 
—-=z(1- “i? + - (sin Ax — Az) My 
3 2! T Jo 








We furthermore abbreviate 


Koo _ K3, = x*y’/2, Kor = K;, + x*/6, 








* ak Qa*y? — 4a°y* . Qr* — 122°y’ 

Nie — 40 ~~ 9 + 9 ’ (42) 
a* — l52*y* — x’y’ 

K > = K ” Neosat -_-—_. 

= . “60 6 


We start the analysis by considering, for the time being, the following problem. 
Find the stress function K of the infinite strip subject to boundary tractions 


*In what follows we omit the subscripts e of A, B, C, D for the sake of simpler notation. The co- 
efficients Ao , By , Co , Do which arise in the K, , J, expressions when n is odd, are given in [2], Eq. (13). 
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5 4 ! 
— -~ z <0 (43) 
0 zr>0 


~ r 7 
or what amounts to the same, find K(z, y) so that 





a, = @, 
44) 
4 9 a2 3.3 dd ( 
2"Ke, = 5 —? | (si “ Me) a. 
2 A 4 =f sin Ax Ax + 31) x° 
It is easy to see that insofar as the first term in the integrand of (44) is concerned 
64 ¢°” si ~ 
2 | Actny + Byshy | ea aD (45) 
wT Jo 


is the solution, A(A), B(A) being given by (40). But (45) is not a satisfactory expression 
since it diverges at \ = 0, moreover no terms are given in (45) which provide the z, 
x*/3! terms of the integrand of (44). We first take care of the singularity of the expression 
(45). Noting that 


1 ’ 


Achy + Byshry = : +rA*(1 - | -sy + —(3 —- y’) 


7 

6! (46) 
yy ~ 2 

_ 6-81 «41 — 66y° + 9y') + oa 


we subtract out the singular part from (45), as shown in the expression (41c). These 
subtracted out terms are nonbiharmonic, so we must add them back in, and this is done 
in the braces { } of K,, . To compensate furthermore for the singularities introduced 
by the terms sin \z/\’ etc. into the braces, we add, under the integral signs of these 
braces, such terms Az, \*z*/3!, etc., that the singularity be removed. The addition of 
z°/6! and z’/2! in the braces outside of the integral sign is finally made in order to 
insure the condition 


abies = 0, |< _ 0 for x > 0. (47) 
For x < 0 there follows then from (41c) 
+ = 0, | _ 2°x*/4! (48) 


Thus, we succeeded in constructing a function, K,, , which satisfies the boundary con- 
ditions (43), (44). Unfortunately, K,, suffers from the defect that it is (because of the 
presence of terms in the braces which are outside of the integral sign) not biharmonic. 
This defect can be corrected by addition of suitable terms, as in (42c), to make it a bi- 
harmonic function K,, ; but then the conditions (47), (48) become (mildly) violated. 

So we start out on a new tack. We consider the problem of finding a stress function 
K,, which gives 


—_ \22°y r<0 (49) 
lo 2«>0 


o* = 0, T 


7Compare Eqs. (11) of [2]. 
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or, what amounts to the same thing, finding a biharmonic function which assumes the 


boundary values 
oS 0, 


2) 2 (50) 
3 *) . 2) ] 
+ 8/42 = | (cos Ax -1¢4 ex) : ‘\ 


(3!) og Jy * 


Proceeding as before, we are first lead to the C and D terms in A,, , and then, noting that 

* ° ‘ ‘ 
Cchrxy + Dyshry , (l1-—y) 4 r (1 y") | “31 3.5 (1 + 3y°) + see (51) 
to the complete expression K,, of (41d) which satisfies the boundary conditions 


lor x > 0: a KP 0 


(52) 
for x <0: Re 2. c=), | ¢ a 82° y/3 


but is nonbiharmonic, and the modified expression K,, of (42d) which does not satisfy 
the condition (52) but is biharmonic. Note, however, that by taking the combination 


a ; K,. (53) 


e ° ye . 6 4 2 ae ° 
we eliminate the nonbiharmonicity of the 6th degree terms 2°, xy", x°y” in the combina- 
tion, and by adding on to (53) a suitable combination of K,, , K., (in the present instance, 
8K.,/9) we eliminate the nonbiharmonicity also of the 4th degree terms. Thus, 


. = ees, e 
K= Ky, - 2 K,, + — Kz, (54) 


9g 


is biharmonic and it satisfies the boundary conditions" 


) > 0 
x, & ’ (55a,b) 
(2%x'/4! r <0, 
) >a _— 
oot: ae. ' (55¢,d) 
—32(x° + x)/9 1 0, 
2 - 
ag a; 0, aa dv° /dy 9 (1 — y)(5 — y’)). (56a,b) 


Our last step is a rather minor one. We take the combination 


63 
= 


K, K — 34K, (K, = K., + Ko) (57) 


as stated earlier, in Eq. (39). In this fashion we modify our tractions to 


( 0 z>0 (58a) 


K? rm => « 
(422 — 13627 <0, 
®The boundary conditions (55b,d) are just as suitable for our purposes as were the proposed con- 

ditions (43). Equations (43) were not an objective, but merely a starting point. The objective is de- 

termination of some biharmonic function K(z, y) which satisfies (55a,c), (56a) and leads to some 4th 
degree polynomial of the type (56b). 
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tof * viet (58b) 
(—1122" + 242 xz <0, 
K?.y = 9, K?..2 = 1 — 8y’ + 7y’*, (58¢,d) 


and thus achieve that the dv°/dy distribution belonging to K, is orthogonal to the 
dv°/dy distribution belonging to K, . This facilitates expansion of a given dv°/dy into 
K° .. functions. Higher K,, functions are constructed similarly. 

We have thus found that the mixed edge value problem of the semi-infinite strip 


o° = 0, dv°/dy = nth degree polynomial in y (59) 


may be converted into the problem of determining suitable distributions o% , r* which 
are 0 for x > 0, and for x < 0 they are n-degree and (n — 1)-degree polynomials, re- 
spectively, when n = even, and (n — 1)-degree and n-degree polynomials, respectively, 
when » = odd. Furthermore, the stress functions of these o% , 7* distributions may 
be derived in a systematic, though very tedious, manner. So one really does not have 
to solve for them, but merely construct their expressions. Nevertheless, for large n 
the Fourier integral representations of K,(z, y), J,(z, y) become unmanageably cum- 
bersome. For the region x > 0 the braces of the K,, , Ku, , Jue , Jnr Expressions [see 
Eqs. (41)] vanish, and what remains may be converted by contour integration into 
rather ‘simple looking” infinite series, as shown in [4], [7], [2]. These are the expressions 
displayed in Eqs. (28) and in the first part of Table I’. 
The mixed end problem of the semi-infinite strip 


7° = 0, u° = nth degree polynomial in y (60) 
leads in a completely similar manner to distributions o% , 7* which are zero for x > 0, 
and for x < 0 they are n — 1, n — 2 degree polynomials when n = even, and n — 2, 
n — 1 degree polynomials when n = odd. The biharmonic eigenfunction expansions 


of the stress functions, corresponding to various distributions u°, are displayed in Eq. 
(29) and in the second part of Table I. 
Once we have the integral expressions of K, , J, we may also obtain, by differen- 


tiation, the boundary values 
% = — Ke ’ o° = Saas . (61) 


Note particularly that these stress boundary values involve only terms from the brackets 
of the type (41) expressions and no terms from the braces. (The terms in the braces are 
required, as shown in [2], for determining the unknown displacement—w°(y) in the 
case of K, dv°(y)/dy in the case of J—along edge x = 0.) By orthonormalizing the 
functions 7° , ¢° into functions ¢, , s, , as outlined in the Introduction, we are finally 
led to the solutions 7, , G, of the “pure’”’ end problem. 


°The minute one tries to use the series in numerical work, the simplicity is replaced by extreme 
complexity. The writer is therefore inclined to believe that the approximate method of “‘self-equi- 
librating polynomials”, developed in [3], and well established in regard to simplicity and adequate 
reliability in [4], [7], [8], [9], [2], [14] is likely to remain the favored technique in the solution of practical 
problems, at least until the time when tables of the functions K,, J, and their derivatives become 
available. (An extension of the self-equilibrating function method to stress problems in polar coordinates 


is outlined in [13].) 
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4. The expansion formula of R. C. T. Smith. While in the previous section we 
carried out the program of solving the mixed and the pure end problems of the semi- 
infinite strip, we did more than was proposed, we obtained, in addition, also the dis- 
tributions o*% , r* which give rise to the boundary values o° , r°. This latter result, while 
of profound mathematical interest is, nevertheless, not germane to the original problem. 
To obtain o% , 7* , i.e., to obtain expressions like Eqs. (39), (41), we had to pay a very 
high price in amount of labor. 


The great amount of labor resulted from the necessity for determining tractions at 


x = —© which balance certain infinite resultants of the o* , r* distributions. These 


contributions at — © were brought into rather sharp focus in [2]'°. However, in our 
immediate objective we are not interested in what happens at x = —o, nor, indeed, 
are we interested in what happens at x < 0. We are interested only in what happens 
to the right of x = 0. It is, therefore, very desirable to find an alternate method, which 
bypasses the determination of o* , r*, and leads directly to the eigenfunction expansions 
of the K, , J, , H, , G, sets. Such a direct route is, indeed, provided in the excellent 
work of R. C. T. Smith [10]. 

Smith in his paper has shown (using a different notation) that the boundary value 


problem of the semi-infinite strip 


4 * — * —_ 

€ = = @& ——— 4 
V'® 0, i 0, (62) 
6°. = given, @°, = given, 


may be solved in the form of a series in the Fadle-Papkovitch functions (27), 


@$=R > C,4,, (63) 
where 
(z, sin® z,/2 cos’ z,) [ (k = even), 
C, = ' (64) 
—(z, cos® z,/2 sin® z,) [ (k = odd), 
Sis . . —¢° 
/ = | [Py + 22.8% , 2,.0%]- ia dy. (65) 
-1 


Os + 20%, 
In particular, for the cases of 


of = 0, ©°.. = dv°/dy = given, (66a) 


7° = 0, ©°, = u° = given, (66b) 


19Ref. [2] gives many details and side lights which could not be compressed into the present paper. 
The author is indebted to Professor Eric Reissner for calling attention to Smith’s work. The papers 
[4], [7], and the first version of [2] were prepared without knowledge of Smith’s paper. The original 
program, and many of the calculations of the present paper were also carried out without the benefit 
of familiarity with Smith’s results. However, the functions J, , Js, Js, Ks, Ks , Ke were obtained on the 


basis of Smith’s expansion formula. 
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formula (65) reduces to 
> +1 +1 
[=a [ ex@r/ay dy, f=-f o2,0? dy, (67a,b) 
. = =i 


respectively. 
For the case (66a), (67a) the function ® of (63) is the stress function (we refer to it 
as function K); for the case (66b), (67b) the z-derivative of ® 


J=@,=@)C,4.. (68) 
is the stress function. Letting 
lm =kll—1)---(l-i+) (69) 
and noting 


] 


a+ +1 
50,,2 se ae ° b+2 


4 cosz , 1 2). 1-2 (4), 1-4 (6), 1-6 
= oe i 21 "2h + 31 "en — 4/ "en T °° (70a) 





— cos’ z,[lz,? — 21% 2,-* + 31 2i-® — -+-J} 

for k, 1 = even, 

o+1 1 n+l 
,0,,1 a: a ° 1+2 
i Oiy dy (1 re 2)(1 + 1) I. a dy 
4 + - - 
oe fee? — 212i? + 31 2p-% — --- (70b) 

2k k 


— sin’ z,[lz,' — 212i? + 3l2i-§ — ---J} 


for k, 1 = odd, we determine the integrals listed in Table V and thereby also the co- 
efficients C, of expansions (63), (68), when dv°/dy, u° are powers of y. 

Note that the expansions obtained will, in general, not converge for x = 0, since 
the conditions of self-equilibration are not, in general, satisfied by the distributions 
y'. However, the expansions are merely the building blocks which make up the complete 
self-equilibrating stress functions of Tables Ia,b; the latter expansions, as Smith has 
shown, do converge. Smith’s procedure therefore permits complete solution of the two 
mixed end problems. In order to go beyond the mixed problem and resolve the pure 
end problem, we have to determine the edge values of the 7° , c° distributions, in the 
form given in Table II. However, even this particular representation can be arrived 
at by an extension of the Smith technique (without prior determination of the dis- 
tributions o* , r* which give rise to it) by merely retracing the steps employed in the 
usual application of the residue theorem. Starting out with, e.g., 





. 4(2 os” sin’ oe 
shui, «af (2 + cos’ z,) sin a (1 rs oot ts) sin zy + y cot % cos ay | 





- z, COS’ 2 k 
TT} > Q(z,) sin zy + @(z.)y cos ZY (71) 
~ k 4 cos’ 2k . 


Q(z.) = -" (3 cos z, + 2 sin 2), B@(z,) = (2 sin z, — 2 COS%), 
k k 
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TABLE V. Integrals of ®py'. (Subscript k of z is omitted for the sake of simpler notation.) The upper 


group pertains to k = even, the lower group tok = odd. The integrals vanish when k + 1 = odd. 
eo ol 
L | ef eoy' dr —*__ I ee wy)? dy 
; y ie (+s 2)( (1 ae 1) k y 7] 
0 | —+4 cos 2° 
2 | (—4 cos 2/2")[z” — 4 — 2 cos’ z] 
| 
4 | (—4 cos z/z°)[z* — 242” + 72 — cos’ 2(42” — 48)] 
6 | (—4 cos z/z*)[e" — 602* + 10802” — 2880 — cos z(62* — 2402” + 2160)] 
8 | (—4 cos z/z'°)[z* — 1122° + 50402* — 80640z* + 201600 


— cos’ 2(82° — 6722* + 201602’ — 161280) ] 


— 
en 
2 
+ 
a 
™M 
n 
— 
ay 
2 
f 
=~ 
x 
at 


122’ — sin’ 2(3z” — 12)] 


| (—4/2" cos z)[z° — 40z* + 3602? — sin’ 2(5z* — 1202’ + 360)] 


ow w 
= 

He 

2 

7) 

N 

| 


2° cosz)[z° — 842° + 2520z* — 201602’ — sin? 2(7z’ — 420z* + 75602? — 20160) ] 


“J 
—~ 
ne 


9 | (—4/2"' cos z)[z'"° — 1442° + 9072z° — 2419202* + 18144002” 
— sin’ 2(9z° — 10082° + 45360z* — 7257602 + 1814400)] 


we may write, on noting 


d 


dz, 


22 2z,) = 2(1 + cos 2z,) = 4 cos x , (72) 


(sin mek ack) 


the integral representation (0 < ¢€ < R(z2)) 





1/2 f° ** @2) sin zy + B(2)y cos z vi 
ogi ae Eke [ ey Y de. (73) 
2at Jesi sin 22 oe 22 
Introducing 
A= -z (74) 


we finally convert (73) to 


} <* A)chry + V(A)yshry Z_ 
es ~ shor oA, (7: 
-? I [= sh2x + 2r + CA, y) | dr, 79) 


where C(A, y) is so determined that the infinities of the U, U terms are compensated. 
This leads us back to our 72 expression in Table II. (This is the way the edge distributions 
0% ,03,03,7%, 7%, 79 Of Table II were determined.) 

It may be added that the approach of Sec. 3, while evidently superfluous in the 
final establishment of Tables I to V, is by no means superfluous in the deeper insight 
it has given for the solution of the end problem. In fact, it was the very search for the 
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distributions o* , r* which give prescribed oa , r° that provided the guiding idea which 

lead to the solution of the pure end problem of the semi-infinite strip’’. 
Acknowledgement. The author is very much indebted to his colleague, Mrs. J. S. 

Born, for the tremendous task of computing the Tables I to V and the (not shown) 

power series expansions of the functions in the Appendix, and for checking all derivations. 
Appendix. The functions A(A), k(A, y), j(A, y)”. 


A, = sh2X + 2X, Ao = sh2X — 2X, 


A.k,, = chdrshrxy — shr-ychry, 


A.k,. = (shy — chr)shyxy — Achr-ychry + 4yA, , 
Aokeo = d*shdchdy — d*chr-yshry — 3(1 — y*)Ao , 
Aokso = (shd — Achd)chry + Ashdr-yshry + 3(1 — 3y*)Ao , 


A.jee = (AchX — shd)chrky — Ashdr-yshry, 

Avjre = (A*shd — 2dchdr)chry — dA*chr-yshry + 4A, , 
Aojeo = (A*shd — d*chd)shdy — A*chr-ychry + ByAo , 
Aojro = (2Ashd — A*chd)shdry + A*shr-ychry — FyAo . 
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or 


“There exists a second type of pure end problem which is of interest, namely, the one where on the 
edge x = 0 both displacements u° and v° are specified. Since this group of problems is usually character- 
ized by infinite stresses in the corners (xz, y) = (0, +1), the present scheme of infinite expansions cannot 
be adopted without first isolating the infinity or taking other precautionary measures, see [11], [12]. 

3The writer shall be glad to pass on, to those interested, power series expansions of the k(A, y), 
j(d, y) functions (these are needed for evaluating the = 0 values of the integrands of Table II) as well 
as tabulated values of the functions for \ = 0, .4, .8, ... , 6.4, aty = 0, .2, .4, .6, .8, .9, .95, 1.0. 
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BOOK REVIEWS 
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Les fonctions orthogonales dans les problemes aux limites de la physique mathematique. 

By Theodore Vogel. Renseignements et Vente au Service des Publications du CNRS, 

45 Rue d’ulm, Paris, 1953. iv + 191 pp. $3.50. 

This is a short compendium on orthogonal functions and their applications in solving boundary 
value problems in mathematical physics. The book is intended for engineers, physicists and applied 
mathematicians. It is not a textbook, but rather an exposition of results and methods. Nevertheless 
the book is very readable, and while most theory is stated without proof, the main ideas are brought 
out clearly and forcefully, and the reader gets a good picture of the subject. An extensive but selective 
bibliography makes the book particularly valuable. An unusual feature of a book directed primarily 
toward applications is the high standard of rigor in stating results. In particular, the Lebesgue inte- 
gration theory is assumed and used. 

The Table of Contents gives an idea of the large amount of information which is contained in this 
slim but substantial volume. 

Chapter I. Orthogonal functions and differential systems. 

1] Closed sequences of orthogonal functions. 

2] Differential systems and orthogonal functions. 

3] Closure of certain sequences of Green functions. 

4] Perturbed systems. 

Chapter II. Study of some particular closed sequences. 

5] Trigonometric functions and Fourier series. 

6] Bessel functions. 

7] Legendre functions. 

8] Orthogonal polynomials numerals. 

9] Mathieu functions. 

Chapter III. Examples of applications. 

10] Separation of variables. 

11] Boundary value problems for the wave equation. 

12] Boundary value problems for the Laplace equation. 

13] The heat equation. 

14] Equations of the fourth order. 


15] Perturbed problems. L. Bers 


Numerical methods. By A. D. Booth. Butterworths Scientific Publications, London, 
1955. vii + 195 pp. $6.00. 

The orientation of the book is best indicated by the following passages from the preface. “Some 
existing books on numerical analysis lay much stress upon the detailed form in which a given procedure 
is to be laid out; --- such concentration upon actual numbers obscures the underlying mathematical 
basis upon which the work rests, and such tabulations are almost entirely absent from this book.”’ 
“The classical methods of hand calculation are, to a greater or less extent, unsuitable for the modern 
machines, and only by having a thorough knowledge of the underlying mathematical principles, is the 
programmer likely to make effective use of the new tools.” The scope of the book is shown by the follow- 
ing chapter headings, each parenthesis containing the number of pages devoted to the topic: The nature 
and purpose of numerical analysis (6)—Tabulations and differences (5)—Interpolation (16)—Numerical 
differentiation and integration (23)—The summation of series (5)—The solution of ordinary differential 
equations (16)—Simultaneous linear equations (38)—Partial differential equations (32)—Non-linear 
algebraic equations (19)—Approximating functions (9)—Fourier synthesis and analysis (10)—Integral 
equations (8). Since the space does not permit anything like an exhaustive treatment of these topics, 
the author restricts himself to discussing a small selection of proved methods. The exposition is clear 
and easy to follow. The illustrative examples frequently have unusual features that challenge the 
W. PRaGER 


numerical analyst. 
(Continued on p. 112) 








THE INTERACTION OF AN ACOUSTIC WAVE AND AN ELASTIC SPHERICAL 
SHELL* 


BY 
P. MANN-NACHBAR 
Missile Systems Division, Lockheed Aircraft Corp. 


Abstract. The effect of the impact of a plane pressure wave on an elastic spherical 
shell is considered, taking into account both the incident and diffracted waves. An 
infinite series (mode) solution is used and numerical results are obtained for the case 
of a steel shell in water. It is shown that a simple quasi-static system can be used to 
find a very good approximation for the stress. The total acceleration on initial impact 
is found exactly and seen to differ markedly from the accelerations associated with the 
lower vibrational modes. 

1. Introduction. This problem was suggested by G. F. Carrier in consequence of 
work previously done by him on a related problem [1]' (see also [2]). 

An attempt had been made to determine the response of a cylindrical shell to an 
acoustic wave. The form of the functions dealt with in the analysis made it difficult 
to obtain accurate explicit results. If the obstacle is taken to be spherical in shape, we 
still have a practical though highly simplified model of an actual physical structure; 
moreover, the problem becomes mathematically simpler, admitting of exact solutions 
for the deformation and accompanying strains in the elastic body. It was therefore 
decided to treat the case of the sphere in detail. 

2. Forced vibrations of a thin spherical shell. Consider a closed shell of thickness 
h with h «< R where R is the radius of the middle surface. The motion of any closed 
oval’ shell, in particular a spherical shell is, by a theorem of Jellett [3], primarily ex- 
tensional. Therefore, the general membrane theory of shells is applicable. If we locate 
the origin of our coordinate system at the center of the shell and choose as the z-axis 
the direction of propagation of the incoming wave, then we have the additional simplifi- 
cation of symmetrical loading (see Fig. 1A). 

The equations of dynamic equilibrium for an element of shell may, therefore, be 
written [4], 

2 


0 ‘ 
= (N, sin 0) — N, cos @ — ph se sin 6 = 0, (2.1) 


ss aw 

Ne+N,- (s — ph owe 0. (2.2) 
Here N, and N, are the normal forces/unit length acting on the sides of the element 
(see Fig. 1B), s is the external load (radial in direction), p is the shell density, and v 
and w are the tangential and radial components of the displacement. 
N, and N, can be eliminated from the Eqs. (2.1) and (2.2) by the use of Hooke’s 
*Received February 10, 1956. The results in this paper appeared previously as Brown University 
report to the Office of Naval Research under Contract N7onr-35810, No. B11-11/34, 1953. 

1Numbers in brackets refer to the Bibliography at the end of the paper. 

*Principal radii of curvature finite and of the same sign. 
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Fic. 1. Geometry of the problem. 


law and the expressions for the strains in spherical coordinates. The resulting equations 
can be further simplified by the introduction of the dimensionless variables (£ is Young’s 
modulus), 


* = 


= w/R; v* = v/R; r = t/t = Et’/(R’p); s* = s/s, = Rs/(Eh). 


w 


After some manipulation, we obtain for the equations of motion the following relations 
in w* and v* only (»v is the Poisson’s ratio). 


M(v*) = 2(1 — v’) ou + (1 —»’*\(1 — ») oe = + Qu* + (1 — ») oe 
—2 oe — (1 — ») ns + 2Qv* cot? 6 + (1 — v)(cot® 6) ~ uk - 
- Kot O + (1 — »)(cot 6) = ov" stil +9] (co! 8) ~ 
a0 00 
— v* esc 6+ a | = (1 -— +) " : (2.3) 
Liw*) = (cot ao + 2x Fj -) ra + od ae 
Liaw , 3 i oe; y= ~ 1 +») a'w* 
OT 
= cot 6 * — (1 — »’) ae te, + (1 — v)s*. (2.4) 


3. Acoustic wave propagation. The linearized theory of wave propagation yields 
the acoustic wave equation 


Ad — d’¢,,. = 0, (3.1) 


where ¢ is the velocity potential, Ad is the Laplacian in spherical coordinates, \* = 
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R?/(tic?) = E/(pc’), and c is the acoustic speed of the fluid. The pressure perturbation 
P is given by 
P 


—p*?¢,,(r, 6, t), 
(3.2) 
p* = p,R*/t, = p,E/p, 


where p; is the fluid density. The applied stress, s, of (2.2) must, of course, be the same 
as the acoustic pressure P at the surface of the sphere. We have, in fact 


s* = —P/s — (p*/so)¢,-(1, 6, T) seed (p,R/ph)¢, (1, 6, 7) 
= B¢,-(1, 6, T). 


4. The interaction problem. We now pose the following problem. An incoming 
plane wave, ¢) , obeying (3.1) impinges on an elastic spherical shell causing it to vibrate 
according to (2.3) and (2.4). The vibrating body acts as a scatterer and, to a lesser 
extent, as a radiator. The outgoing waves (scattered and radiated) also obeying (3.1), 
in turn influence the nature of the vibrations. At the surface of the shell, the radial 
velocity of the fluid, ¢, r(1, 6, 7) must be equal to the radial velocity of the shell, w*,(6, 7). 
We wish to determine the motion of the sphere and the pressure distribution associated 
with the incoming and outgoing waves. 

The pressure associated with the incident wave is taken to be [1] 


(2 exp [62 — r/A)] 2S 7/A 


Il 


P= (4.1) 
lo z>1/r 
so that the initial velocity potential is 
(Qod /p*5) exp [5(z — 1/d)] — (QodA/p*5d) z2<1/xr 
b> = (4.2) 
0 z> 1r/x. 
For simplicity let Yo , ¥, x, W and V be defined as: 
do = (Qo/p*) vo , & = (Q./p*)(vo + Y) = (Qo/p*)x, 
(4.3) 
w* = (Q./p*)W, v* = (Q./p*)V. 


Then our boundary value problem is defined by the following set of equations’. 
L(W) = Blx.-(1, 0, 7) cot 6 — (1 — v*)x,srr(1, 8, 7) 
+ X.eor(1, 6, 7) + (1 zai v)x,-(1, 6, 1), (4.4) 


M(V) = Atl - v*)x,er(1, 6, 7), (4.5) 

Ay ta so = 0, (4.6) 

W..(@, 7) = x.-(1, 8, 7), [boundary condition] (4.7) 
(\/6) exp [8(z — 1r/d)] — A/5 z<1/d 

Yo = (4.8) 
0 z>1/xr. 





*The operators L and M which appear here are those given in (2.3) and (2.4). 
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These equations will be more easily handled if Fourier transforms are first introduced 
to eliminate the time dependence. Denote the transform of a function F by F’. Then 
F and F’ are related by 

yo e ’ . 
F'(r, 0,é) = | F(r, 6, r) exp (—iér) dr, 
hs (4.9a) 


Imé& = —1a, a> 0, 


Fir, 0,7) =: 


‘ 


= / F'(r, 0, &) exp (itr) dé. (4.9b) 


bo] 


Applying (4.9a) to (4.4) through (4.8) we obtain 
L"(W") = iBe[x", cot 6 + x70(1, 0,€) + (1 — x” + (1 — »*)x’], (4.10) 


where 


a 2 — 


Ps (9 _ 2 >) | ae a — 
L (cot a(2 E56) +2+ ag ~ § 29? (1 v’)é 


— (1 — ve? + 2(1 + ve’, 


M"(V") = iBt(1 — v*)x” , (4.11) 
where 
M’ = -201 — 2 + (1 — Y)(1 — vet + Q& — (1 — VP 
9” ‘ , <u 
— 2—5 + (1 — ve + 2 cot? 6 — (1 — ve’ cot’ 0 
06 
— 2cot 8) Ld + (1 — v)é*(cot 6) <= +(1+ | ot 8) a 6+ *I 
C08) 80 ve {cot 0) ae [ — = ae 
Ay’ + r\{,#y’" = 0, (4.12) 
i¢W"(0,&) = x(1, 0,8), (4.13) 
, _exp (—12N€) = a ae eo 
y iEI(8/) + 18] f() exp (—aré cos 8); 


SE) = {2&[(6/d) + 2€]}. (4.14) 
It can be shown, using the method of separation of variables on (4.12), that y” is 
of the form 
rLJ,.4(dér) + CY,.,(dér)]P,(cos 8). 


Since the waves associated with ¥” must be outwardly moving, take C = —i and write: 
v7 = >> ATER (rEr)P,(cos 8), (4.15) 
where 


hy? (rér) = (9 /2dér)* [Jana y(Abr) — 2¥ 44 (AEN). 
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y> may be expanded similarly 


Vr, 0,8) = —f@ ¥ (n + 1)(—9"j,08)P,(cos 6), 


n=0 


(4.16) 
= > B,(t)j,(Atr)P,(cos 8), 


n=0 
where 
ja(dér) = (4 /2dEr)*J,,.\(AEr). 
From the work of Lamb [5], we know that the complete solutions for W* and V7 may 


be written: 


W'(0,t) = >> W7)P,(cos 6), (4.17) 
V7(6,8) = — > Vi®P\(cos 0) = — > V7P%(cos 6) sin 6. (4.18) 


Substituting these expressions into (4.10), (4.11) and (4.13), we get three algebraic 
equations for the three unknowns W? , V7 and A? , which are readily solved to give: 











r______—(2n + _1)(—7)"C, 7 
aiken AEE) [(6/A) + AEC AR?” + iCALY]’ (4.19) 
yr —(2n + 1)(—1)"(1 + »)8 use 


"  NEG(EA) + UCM” CA”)? 


Crh’ + iC, ho” 
F T —— 1 J 1 J 
a i T CMe” + sh 





|ke + 1)(-7)", (4.21) 


where 

C,=1-)@’4+n)-24+0-Y)24+ 0 — vd! — 21 + ve, 

C, = ipl — YP — (rn? +n) + (1 — J. 
The quantities of physical interest are the stresses, radial acceleration and the total 
pressure distribution. These can now all be found, at least in principle, from (4.19), 
(4.20), (4.21) and the inversion formula (4.9b). 


5. Numerical results for the shell. For each vibrational mode, the transforms of 
the stress components will be linear combinations of W? and V2 ; e.g. for the nth mode, 


N% will be given by 
1 
[(Q.R/(1 — “all « + »)WP,(cos 6) — V™P\(cos 6) cot 6 — vV7 2Paloos 0) ) 
The transform of the total radial acceleration, d°w/dr’, will be —(Q.R?/Ehg)?W’" or, 
for the nth mode only 
—(Q)R?/Ehp)#?W’.P,,(cos 8). 


From (4.19) and (4.20) we see that these expressions are regular in & except for a finite 
number of poles. The theorem of residues and Jordan’s lemma can therefore be used to 
evaluate the integral of (4.9b). The computation was carried out in detail for the three 
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lowest modes, and the results are shown for the point 6 = z in Figs. 2 and 3. N, only 
is shown in Fig. 2 since the two tangential stress components are equal for 06 = rm. 
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Fic. 2. Hoop stress at @ = = for the first 3 modes. 


The parameters were taken as 


B=17, 6=0, W=25/3, v=3, 


the values appropriate to a step wave of infinite length impinging on a steel shell in 
water. 

From Fig. 2 it becomes evident that the series representing the stresses in the shell 
converge rapidly enough so that a very good approximation to the total stress may be 
obtained by considering just the first few modes. Of these, the zeroth mode, which 
has the steady state solution as its limiting form, is the most important. In this con- 
nection it should be noted that the curves in Fig. 2, which were made for 6 = x, show 
the largest stresses which can occur in the first and second modes at any time. The 
zeroth mode is, of course, @ — independent. 

The computation of the acceleration is not so simple. A picture of the total accelera- 
tion cannot be gained by looking at the lower modes. In fact, it would seem from Fig. 
3 that a few terms of the series }>*., 0°w,/dr’ P, (cos @) will not approximate with 
acceptable accuracy the total acceleration for all r and @. The rate of change of radial 
acceleration is very great at 7/A = —1, 6 = a; probably the total radial acceleration 
will be discontinuous at r/A = —1, 6 = z. That this is actually so can be seen by re- 
calling that the incoming wave reaches the point (1, 6) on the sphere at time r/A = cos 8. 
In particular, the point of initial impact, 6 = 7, is reaci.od at r/A = —1. This means 
that the pressure and, therefore, the radial acceleration as well are discontinuous and 
that the series representing them do not converge uniformly for r/A = cos @. The solutions 


1957) INTERACTION OF ACOUSTIC WAVE AND SPHERICAL SHELL 


























89 
DYNAMIC CASE 
---—- — QUASI-STATIC CASE 
1.0} Fan 
i \ 
= J \ 
/ 
4 \ 
ery , 
1 
Ra, f 1 
\ 
Hi! ‘ 
/ \ 
1 
> ’ \ 
| t 
iN 
3 \ 
' a H 
if ' 
Oo 10 \20 _,30 
T/A+1 a/aal 
-tF 
-2- -2b -2+ 
aa -3F 
7 ZEROTH MODE » FIRST MODE “4h 





SECOND MODE 
Fic. 3. Radial acceleration at @ = x for the first 3 modes. 


for these points must be found in closed form, i.e. when the transform, ¢?W” (4.17), is 
inverted, the order of summation and integration cannot be interchanged. This difficulty, 
associated with the use of the series expansion as a method of solution, will be encount- 
ered again in Sec. 6 when the resultant pressure distribution is discussed. By anticipating 
the results of that section, we can find the initial value of 0’w/dr’ at 6 = =. 


Equation (2.2) gives, for r/A = —1 


h a°w Eh aw 
— p a mma | 


at? sR? ar 
—4s —4Eh dw 


~” 


NI 


’ 





ae Grew 
Taking 


s = —2Q,[(6.6), 6 = x] 


we have 


- 4Eh dw 


OR a 





Thus, while we can obtain a very good approximation to the stresses by considering 
only the first few modes or even the lowest mode by itself, the same is definitely not 
true of the acceleration. The actual initial acceleration is about 40 times as great as the 


maximum acceleration in the zeroth mode. The agreement obtained by considering the 
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first and second modes along with the zeroth is not appreciably better; the results still 
differ by more than a factor of 6. 

Quasi-Static Case. It is of interest to compare the hoop stresses for the zeroth mode 
with those for the quasi-static case; i.e. for the case in which both scattering and the 
inertia forces due to deformation are neglected. 

Since the effects of only the incident wave are considered, the pressure at time ¢ 
is given by 

QodaR(te + R) _ Qol(z/A) + 1) wit thes 
4nrR 2 





pressure = 
Qo 1 < 7r/X. 


The stresses are given by 





j=Peule» SN ose ae. 


—2N,/(RQ.) appears as the dotted line in Fig. 2, where it can be seen that the very 
simple quasi-static case approximates the more exact dynamic case very closely. The 
stress developments, while slightly out of phase, are essentially parallel with a difference 
in maxima of only 10 per cent. 

For —1 < r/A < 1, there will be an unbalanced force, due to the incident wave, 
acting on the sphere. This will result in a rigid-body acceleration in the z direction. 


The magnitude of the force is 
Q.rR*[1 — (r/d)*] 
therefore 


4rR’*hp os = QorR*[1 —(r )*] 
OT 


or 


4Eh d°2 <2 
QoR? az? = | — G/N 
This has been plotted for purposes of comparison in Fig. 3. 

The formal resemblance is greatest for the first mode. The maximum acceleration 
in the zeroth mode is only one fifth the rigid-body acceleration. For higher modes, the 
maxima move closer to the rigid-body value of 1. By referring to the previous section, 
however, we see that the maximum total acceleration associated with deformation is 
8 times that of the rigid body. 

6. Resultant pressure distribution. From (3.2), (4.1), (4.3), (4.15), and (4.21) 
and (4.9b) the total pressure is known to be 
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Protas = P; (incident) + Py, = —p*d,, = —Qx.- ; 


ED nae ’ T 
= Q— Se [exp (env 'e ae, 





Il 


Oe es os te ; 
05 + fg exp (igs) SK MAR+ D jmragn| 


n=0 





C.aAor + iC ong 
aT ag ORE P,(cos 6) dé. (6.1) 
1 n 2! 'n 


Unfortunately, the order of summation and integration in (6.1) cannot be inter- 
changed for all @, r and r; the total pressure cannot be approximated by the sum of the 
first n vibrational modes. This can be seen for the specific case r = 1, r = —A, 0 = 
as follows. At the moment of initial impact, 7 = —\, we have for each mode and for all 
0, N, = N, = 0°v,/d7*? = 0’w,/dr* = 0. The equations of equilibrium, (2.1) and (2.2), 
require that the pressure at the surface of the sphere must likewise vanish for each 
mode at r = —\A, so that the total pressure would be zero for all 6. We should expect, 
however, from what is known of the theory of scattering of plane waves, that the pressure 
(, would be doubled and not reduced to zero. 

vy” must be found in closed form if P7,,,; is to be evaluated. It has not proved feasible 
to carry out the summation for all 6, r and r. However, it was noted that the expansions 
for ¥5 and y’” are very similar for r = 1, and & very large or [(r/A) — cos 6] very small. 
This fact can be used to obtain the pressure at the surface of the sphere for r/A = cos 0. 

By referring to the expression for ¥5(4.16), it is seen that if we set AE = ¢, the pressure 


associated with the incident wave may be written 


P,=Q=- Gos [- Sty > j(¢r)P,(cos 6)(2n + 1)(—1)" df. 
J —tar-o@ $ n=0 


For ¢ > Band r = 1, Py,, (see 6.1), becomes, to first order 


Qo fi" exp (ite/N) 


ee a § n=0 





_ sin [¢ — 2 — 2/2] (on + 1)(—i)"P,(cos 0) dt (6.2) 


and P, reduces to 


— Voi “F@4*2 oxy (itr/d) = 


eee c n=0 


P, = 





_ cos [¢ — orl) — 1/2] (2n + 1)(—1)"P,(cos 6) df. (6.3) 


Equations (6.3) and (4.14) tell us that 


t exp (—7f¢ cos @) > > cos [¢ — (nw/2) — 2/2](2n + 1)(—1)"P,(cos 6) 


f[-o n=0 
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so that to first order: 


Qe f mr exp {i¢[(7/A) — cos 6]} 


Pu=35 exp [t(1/2) cos 6] df, 
oT J-jah-@ $ 
opi J (6.4) 
Q, exp [i(a/2)(cos @ + 1)] for 0 < [(7/A) — cos 6] « .05, 
= 0 for [(r/A) — cos 6] < 0. 
Higher order terms may be obtained in the same way. To second order: 
Q, f-'™*® exp fitl(r/r) — cos 6]} f 7 
Py = = [ o=P isl(r x — L exp [a(x 2) cos 6] + Bat 
amy. fark=—« § t 
+ i exp (i(7/2) cos 6)] — 5: exp [i(7/2) cos ah dt, 
<5 
= Q,{exp [7(r/2) cos 6 + ix/2] — [(r/A) — cos 6][1 + B][I1 (6.5) 


+ exp (i(7/2) cos 6 + in/2)] + 2/2[(r/d) — cos 6][exp (i(r/2) cos 6)]} 
for 0 < [(7/A) — cos 6] < .05, 
= 0 for [(r/A) — cos 6] < 0. 


Additional terms will be of little value since the representation is valid only for 
¢ > 8B, [(r/A) — cos 6] «< .05. The terms found so far, however, are sufficient to tell us 
some things of importance. 

The incoming wave will reach the point (1, @) on the sphere at time r/A = cos 0. 
The pressure on impact for each 6 is given by* 


Proar = P; + Re Pr, = Qo + Qo cos [(x/2)(cos 0 + 1)). (6.6) 


At @ = 7m, the outermost point of the sphere, the first effect is that of a plane wave 
hitting a rigid wall and we have 


Py = P, ; Protar = 2Qo - 
At 6 = 2/2, the wave just grazes the sphere and therefore 
Pi - 0; Protai = Qo ° 


As 6 varies from z to 7/2, the initial pressure varies continuously from 2Q, to Qo . 
If the sphere were rigid, the steady-state pressure distribution would be given by 


Prrota - P; 0 < 6 < T. 


The results of Sec. 5 indicate that we will have asymptotic values RQ,/2 for the 
stresses and zero for the radial acceleration, which also correspond to a uniform pressure 


of Q = P;. 





‘The elastic waves in the shell will travel more rapidly than the acoustic wave and will result in a 
pressure, P = 0, at (1, 0) before the time r/A = cos 0, 6 ~ x. However, this effect is negligibly small 
compared with the one we are considering. 
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At present, this is about all that can be said on the subject of the pressure distri- 
bution. While the series method is very convenient for finding the stresses and may be 
adapted to give results for the radial acceleration, it does not lend itself to an analysis 
of the pressure distribution. It would seem that, for this aspect of the problem, another 
approach is called for. 

7. Remarks. In a paper which appeared after this work was completed, J. H. Huth 
and J. D. Cole [6] consider the related problem of the effect of a shock wave on an air- 
borne elastic spherical shell. 

This problem is more difficult than the one treated here, and the authors have 
introduced the simplifying assumption that the effect of the diffracted wave on the 
shell may be neglected, or, in other words, that the effect on the applied pressure of 
both the motion of the sphere and the sphere itself may be neglected. 

The stresses, but not the radial acceleration, are computed. The problem of finding 
the resultant pressure distribution does not, of course, arise. 

It would be of interest to compare their results for the limiting case of a very weak 
shock (i.e. speed of shock wave — speed of sound) with the stresses that would be 
obtained here by using for 8 and \ (see Sec. 3) the values appropriate to a steel shell in 
air. Such computations are now under way. 
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—NOTES— 


ON THE CONDITIONS OF VALIDITY OF RIEMANN’S METHOD OF 
INTEGRATION* 


By AUREL WINTNER (The Johns Hopkins University) 


1. Introduction. The traditional treatment of Riemann’s classical formula depends 
on certain additional assumptions which are not necessary for the existence of a solution 
(Sec. 2) but, being based on the adjoint equation (Sec. 3), are necessary for the customary 
introduction of Riemann’s Green function. 

The purpose of this note is to point out the resulting complication and to show that 
it can be overcome by the application of a simple device in the proof. The main result 
is that italicized in Sec. 8. 

It will be sufficient to consider only one of the standard cases, the case in which the 
boundary values of the unknown function itself are assigned along a path consisting of 
two characteristics which meet at a point. In fact, it will be clear from the nature of the 
arguments to be applied that all considerations remain valid for the case in which the 
data are, for instance, Cauchy data proper, the case in which the unknown function 
and its normal derivative are assigned along a continuously differentiable path on which 
no direction is a characteristic direction. 

2. Picard's theorem. On the closed rectangle 0 S x S ¢,0 S y S 9, which will 


be denoted by Q = Q(é, n), consider Riemann’s hyperbolic partial differential equation 


Us, + a(x, yu, + W(x, yu, + c(z, yu = 0, (1) 
with the boundary data 
u(x, n) = ¢(2), u(t, y) = Wy), (2) 
where ¢(z) (0 S x S £) and Y(y) (0 S y S 7) are given functions satisfying 
ef) = Wn). (3) 


The appropriate conditions, to be imposed on (a, b, c) and (gy, ), respectively, are 
as follows: (i) the coefficient functions a(z, y), b(z, y), c(z, y) are continuous on Q and 
(ii) the boundary functions g(x), ¥(y) are continuously differentiable. In fact, the situa- 


tion is as follows. 
The pair conditions (i)-(ii) assures that there exists on the rectangle Q a unique 
function u = u(z, y) having the following properties: 


u, ,u, and u,,(= u,,) exist and are continuous, (4) 


the relation (1) is an identity on Q and the two equations (2) are identities on the sides 


y=12=£ofQ. 
This will be referred to as Picard’s theorem. It can be proved by the method of 


*Received February 3, 1956. 
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successive approximations and is the precise formulation of what is actually proved in 
Picard’s writings on the subject.’ 

3. The formal adjoint. In contrast to what is supplied by Picard’s theorem, Rie- 
mann’s method of integration replaces Eq. (1) by its adjoint, 


Vey — [a(z, y)v]. — [d(z, y)v], + c(z, y)v = 0, (1 bis) 
to which certain boundary data 
v(x, n) = p(x), —v(é, y) = oly) (2 bis) 
satisfying 
p(—) = a(n) (3 bis) 


are assigned. 

The traditional presentation of Riemann’s method’ depends on an application of 
Picard’s theorem to the case in which the Eq. (1) and the boundary condition (2) are 
replaced by the adjoint equation (1 bis) and a certain boundary condition (2 bis), re- 
spectively. A moment’s reflection shows, however, that this application is inadmissible 
under the assumptions (i)-(ii) placed on (1)-(2) by Picard’s theorem. For, since (i) 
merely assumes the continuity of a(z, y) and b(z, y), the derivatives a,(z, y), b,(z, y) 
need not exist (or, if they do, they need not be continuous) and so Eq. (1 bis) cannot 
in general be expanded into 


Vy + A(z, y)v, + B(x, yv, + C(z, yv = 0 


with certain functions A, B, C, still less with continuous functions A, B, C, as it is re- 
quired by (i) when Eq. (1 bis) is identified with Eq. (1). 

4. Methodical remarks. The purpose of this note is to point out an easy way out 
of the implications of this predicament. 

For the case of the elliptic analogues of the hyperbolic equation (1), a similar difficulty 
was circumvented by Lichtenstein,*® by using an “integrated form” of the adjoint. In 
the hyperbolic case at hand, a simple and direct approach, which makes explicit a 
function space allowable for the adjoint under the original assumption (i) on Eq. (1), 
can be obtained. 

Today, Lichtenstein’s result can be interpreted as a manifestation of L. Schwartz’s 
distribution theory*. From the point of view of the theory of distributions, the success 
of the explicit approach in the hyperbolic case will center around the following formal 
fact. If the set-function I*(f) (of the set q) is the integral of a continuous f = f(z, y) 
over a rectangle gq = (x, Sz S 2%,Y; S y S yz), and if f is one of the (continuous) 


See, e.g., E. Picard, Legons sur quelques types simples d’équations aux dérivées partielles avec des 
applications 4 la physique mathématique, 1927, pp. 123-134. For an extension to the case of certain non- 
linear equations, a case the treatment of which is not based on the method of successive approximations 
but on that of equicontinuous functions, see P. Hartman and A. Wintner, Amer. J. of Math. 74, 836-843 
(1952) 

*See, e.g. E. Picard, op. cit., pp. 147-151 or G. Darboux, Legons sur la théorie générale des surfaces, 
2, 71-81 (1889). The literature consulted on the extension of the classical conditions on Riemann’s 
method (see, e.g., G. S. S. Ludford, J. of Ratl. Mech. and Anal. 3, 77-88 (1954), where further references 
are given) does not go into the problems of the adjoint considered in this paper 

*L. Lichtenstein, Bull. Acad. Polon. des Sciences 1931, 571-598 

‘L. Schwartz, 7'héorie des distributions, vol. 1-2, 1950-1951 
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derivatives z, , z, , 2, of a function z = 2(z, y), then the evaluation of J*(f) involves 
no differentiation of z(z, y) in any of the three cases f = z, , 2, , 22, (it involves integrations 
in the first two of the three cases). 

5. Reformulation of the adjoint problem. Suppose that the coefficient functions 
a, b, c of Eq. (1) are subject only to the continuity assumption (i) on Q. Then Eq. (1 bis), 
as it stands, is meaningless in general. But it appears in a meaningful form if it is inte- 
grated over the rectangle having (x, y) and (é, 7) as opposite vertices, where (zx, y) is 
any point of the given rectangle Q = (0 S x S$ ¢,0 S y S 7). In fact, the formal result 


of this integration is 


v(é, n) — v(€, y) — v(x, n) + v(2, y) 


i nf 
dt — | b(s, t)v(s, t) | 


hice 


~ s= 


~ [ a(s, t)v(s, t) 


vy s=z 


ds + [ / c(s, t)v(s, t) ds dt = 0. 


In view of the boundary data (2 bis) and the condition (3 bis), this can be written in the 


form 
en 
c(s, t)o(s, t) ds dt 


v(x, y) = —p(€) + p(x) + ofy) - [ | 
Py pé , 
+ [ {a(é, t)v(é, 1) — ala, dv(ax, t)} dt + | {b(s, n)v(s, n) — O(s, y)v(s, y)} ds, 


[where — p(é) can be replaced by — o(n)]. Let the latter formula for v(z, y), a formula 
which is an integral equation for the unknown 2, with a, b, ¢ and p, o as data, be referred 
to as Eq. (1*). 

Under appropriate assumptions of differentiability, Eq. (1*) is equivalent to Eq. 
(1 bis) and (2 bis) together. But Eq. (1*) is meaningful under the following pair of 
assumptions also: (i*) the continuity assumption (i) of Sec. 2 for a(z, y), b(z, y) c(z, y) 
on @ and (ii*) the continuity of p(x), o(y) on the respective intervalsO S25 t,08 
y = 7 [note that (ii*) requires of p, o less than (ii) in Sec. 2 requires of ¢, y]. 

6. The dual of Picard’s theorem. Corresponding to the circumstance that the 
formulation (1*) of the adjoint problem is free of any differentiation, it is natural to extend 
the solution class, defined by condition (4) in Picard’s theorem, as follows: a “solution” 
v = v(z, y) should mean any function defined on Q in such a way that 


v(x, y) is continuous (4*) 


on Q and satisfies the Eq. (1*) as an identity on Q. The fundamental fact, representing 
the true dual of Picard’s theorem, can then be formulated as the following theorem. 

If (1) and (2) are replaced by (1*) and, at the same time, (i), (ii) and (4) are replaced 
by (i*), (ii*) and (4*), respectively [which means that (i) is retained but both (ii) and (4) 
are extended to the assumption of mere continuity], then there exists on Q exactly one 
solution v(z, y). 

The proof of this dual of Picard’s theorem depends on successive approximations. 
To this end, let »(z, y) = 0 on Q and, if v,(z, y) has already been defined on Q asa 
continuous function, let v,,,(z, y) denote the function which results if v is replaced by 
v, in the sum which represents the expression on the right of Eq. (4*). It then follows 
that v,(x, y), v.(z, y), --* is a sequence of continuous functions which tend to a limit 
function uniformly on Q and that, if v(z, y) denotes the limit function, then v(z, y) is a 
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solution of Eq. (4*). The details are of a routine nature and will therefore be omitted. 
The assumption that there exist two distinct solutions v(z, y) leads, again by successive 
approximations, to a contradiction in the usual way. 

7. The reciprocity relation of the Green functions. Under assumption (i) [which is 
identical with assumption (i*)] choose ¢, ¥ and p, o as follows: 


E ” 
o(z) = exp [| —0(s, 9) ds, vy) = ep [ —ale, 0 at (5) 


and 


E ” 

i sen / b(e, 9) ds, aly) = exp / a(t, t) dt (5 bis) 
(so that p = 1/g, o = 1/). Then conditions (ii), (3) and (ii*), (3 bis) are satisfied. Hence 
both of the above theorems (those of Sec. 2 and Sec. 5) are applicable if a, b, c are just 
continuous. With reference to the rectangle Q, having (0, 0) and (£, 7) as opposite vertices, 
let G(x, y) = G(a, y; £, n) denote that solution u(z, y), subject to condition (4), of Eq. 
(1) which belongs to the Cauchy data (5), and let H(z, y) = H(z, y; &, 9) denote that 
solution v(x, y), subject to condition (4*), of Eq. (1*) which belongs to the Cauchy data 
(5 bis). Then, if (x’, y’), (2, y’’) is any pair of points contained in the rectangle Q, 
Riemann’s reciprocity theorem 


H(z’, y';2",y") = Ga" "_,y";2',y’') (6) 


is valid. 

Since a, b (and c) are now just continuous, the classical proof of the symmetry re- 
lation’, a proof based on Green’s formula, fails to apply. One can however conclude as 
follows. 

It is readily seen that the successive approximations leading to u = G(z’, y’; 2”, y”’) 
and v = H(z’, y'; z’’, y’’) are uniform in the points (z’, y’), (x”’, y’’) of Q and in the index 
k(= 1, 2, ---) together, if (a, , b, , c,), (2 , bo , C2), «-+ , (Qe, Oe , Ge), - ++ Are Sequences 
of functions which are continuous on Q and tend to (a, b, c) uniformly on Q. Hence, the 
relation (6) is true for (a, b, c) if it is true for every (a, , b , c,). Since it is true in the 
classical case, it follows for the case of just continuous coefficient functions (a, b, c); in 
fact, a,(z, y), b,(x, y), c.(x, y) can be chosen to be polynomials. 

8. The extended validity of Riemann’s formula. Without any reference to the 
relation (6) (which will not be used directly), consider again the problem of Sec. 2, 
represented by Eqs. (1), (2) and (3), with (i) and (ii) as assumptions. According to 
Picard’s theorem, this problem has on Q a unique solution, 


u = u(r, y; a, b,c, ¢, W) (7) 


satisfying condition (4). Riemann’s classical formula represents this unique solution 
(7) in terms of the Green function H(z, y; §, ) of the adjoint equation. But the classical 
definition of H is meaningless under the present assumptions (see Sec. 3). On the other 
hand, if the continuous function H(z, y; £, n) is defined, not in the traditional manner, 
but in the way specified in Sec. 7 (as supplied by the existence and uniqueness theorem 


See, e.g., G. Darboux, op. cit., p. 81 
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of Sec. 6), then Riemann’s integral representation of the solution (7) remains valid under 
the assumptions, (i) and (ii), of Picard’s theorem alone. 

In fact, Riemann’s integral representation’ of the solution (7) contains only the 
following elements: (a@) the boundary data g(x), ¥(y) [which are subject to condition 
(3)] and their first derivatives and (8) the function H(z, y; é, n) and their first derivatives 
on the boundary; cf. (5), (5 bis). But item (8) does not involve the second derivatives 
of the function H (derivatives which do in general exist). On the other hand, not only 
the functions g(r), ¥(y) but also their first derivatives, introduced by item (a), are 
controlled by assumption (ii). Hence it is clear that the italicized assertion, concerning 
the general validity of Riemann’s formula, can be concluded by the same argument 
(this time by approximating a(z, y), b(z, y) and ¢(z), ¥(y) as well as de(x)/dx, dy(y)/dy 
by sequences of smooth functions) which was applied in Sec. 7. 


NOTE ON THE AERODYNAMIC HEATING OF AN OSCILLATING 
INSULATED SURFACE* 


By STEPHEN H. MASLEN anp SIMON OSTRACH 
(Lewis Flight Propulsion Laboratory, National Advisory Committee for Aeronautics, Cleveland, Ohio) 


The effect of disturbing the thermal equilibrium of an oscillating conducting surface 
and its surroundings by changing the isothermal surface temperature at a given time 
was investigated in Ref. [1]. It was shown therein that the heat transfer and the thermal 
state of the fluid associated with the oscillating surface can be significantly different 
from that for conduction from a stationary surface with the same initial temperature 
difference. To complete this study it is appropriate to investigate the effect of insulating 
the surface at a given time on the equilibrium state. 

Accordingly, consideration is given herein to a doubly infinite plane surface which is 
oscillating axially (i.e., longitudinally) in a viscous and heat-conducting fluid. It is 
assumed that sufficient time has elapsed so that an equilibrium state exists in which a 
periodic motion of the fluid has been established, and the heat obtained by viscous 
dissipation is all conducted through the surface so that the temperature does not in- 


crease indefinitely with time. In this state the fluid velocity is given by [2] 


u(y, t) = U exp [—(n/2v)'”’y] cos [nt — (n/2vr)'*y] (1) 


and the temperature is [1] 





U*Pr . 1/2 ] : 1/2 y_\1/2 
T.=T.- ry {exp [—(2n/v) **y] — —— [exp {—(n/a) ““y} cos {2nt — (n/a)’’*y} 


— exp {—(2n/v)'”*y} cos {2nt — n/s)'*uihh, (2) 


where wu is the fluid velocity component parallel to the surface, y is the coordinate normal 
to the surface, ¢ denotes time, U and n are the amplitude and frequency of the surface 





§See, e.g., G. Darboux, op. cit., formula (16) on p. 80 
*Received February 20, 1956. 
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oscillations, v is the kinematic viscosity coefficient, Pr is the Prandtl number, c, is the 
specific heat at constant pressure, 7 denotes temperature, a is the thermal diffusivity, 
and the subscripts e and © denote conditions at equilibrium and far from the surface, 
respectively. 

From Eq. (2) it can be seen that the thermal boundary condition at the surface 
(y = 0) for equilibrium is 





U*Pr (3) 


T(0,)=T.. = T.- dc, 


In the problem treated herein the equilibrium is to be disturbed by replacing, at time 
it, , the isothermal condition specified by Eq. (3) by one for an insulated surface 


TO, 1) = 0 t(>&., (4) 

where the subscript denotes partial differentiation. The other boundary conditions are 
T(~,t) = T. (5) 

T(y, to) = Ty, b) (6) 


and the corresponding differential equation is: 


T,.-cdl.= ~ (u,)”. (7) 


Using Laplace transforms, the solution of the problem defined by Eqs. (4) to (7) is found 


to be: 





72 1/2 
rly, ) = Tay, 0) + OE (BY ty erfely/2att — 11") 


+ 2[a(t — t)/x]'” exp [—y*/4a(t — t)]} 


U’ (n /2n)'”? . -1/2 
4: 2, 1 + (2/Pn™ [ cos (2nr + 2/4) exp [—y’/4a(t — 7](t — 1) dr. (8) 
The profiles associated with the oscillations are shown in Fig. 1 for two special cases for 
air. The equilibrium profile (I) is one in which the temperature monotonically increases 

















' 
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so 
T-Ty 
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Fic. 1. Temperature profiles. 
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from the wall to free stream. When the wall is insulated, the wall temperature rises 
[roughly as (¢ — 4))'’*] giving a profile (II) which has a temperature minimum in the 
boundary layer. For larger times the temperature minimum will decrease in magnitude 
and, presumably move toward the outer edge of the boundary layer. Further, as can be 
seen from Eq. (8), the thickness of the temperature boundary layer increases linearly 
with (¢ — t)). 

Although this solution was developed for an incompressible viscous fluid with con- 
stant property values, it is equally valid in the case of a compressible viscous fluid if 
the boundary-layer assumptions are made, if the Prandtl number and the product of 
the density, p, and absolute viscosity coefficients are assumed constant and if y is re- 


placed by 7 where 


ay 
p 
9 = | — dt. 
70 Po 


Under these assumptions the compressible boundary-layer equations reduce to those 


for an incompressible fluid. 
The most important property of an insulated surface is perhaps the recovery factor. 


Using Eqs. (3) and (8), this is 








_170,t)-—T._ _~Pr _ / 1/2 
= Uk, = . + [2Prn(t — t)/r] 


/ 





| [2n(t—to)]*/2 : 
+ i + @/P)™) I cos (2nt + /4 — 6°) dB. 


Since the time dependence itself is not of primary importance, the recovery factor is 
averaged over a cycle to yield 
k+n/ 


mf rat = S$ 2 (PH I[1 + nk — &)/x]”? — Inlk — t)/n}} 


r=- 


Ww Jk 


__ sin (2nt, + m/4) _ F ba foy/2 __ on ;_y1/2 
r(2)7[1 + (2/Pr)'”] {{1 + n(k to) /r] [n(k to) 1] } 
1 [2n(k—to)+24]3/2 
ae creer aed 31 C /4 — B’) dB. 
+ SE OP esacunnn 82 Onk + 2/4 — 6°) ds 


Comparing this result with that for a semi-infinite plate in steady flow (for which r, = 


Pr'’*) and writing the integral in terms of Fresnel integrals yields 


1 yn , (2) 3/2 3/2 
=— 5 (Pr) + 5 {{1 + n(k — t)/e]~° — [n(k — t)/e]*} 





|~ 


_1=-= 


“ 
. 


ee 

Oa[l + (Pr/2)2] (Lh + mk — bo)/e] — [ne — to)/m]"} 
2[l+n(k—to)/w]*/* - 
cos 5 A dg 


2(n(k—-to)/a]t/ 


1 a? / 
+ 4n(1 + (Pr/2)*] (sin (Qnk + 2/4) 


2[1+n(k—to)/]*/2 


— cos (2nk + 2/4) [ sin 1B’/2 is}. (9) 


#2(n(k—to)/e)*/2 
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Obviously, the recovery factor with oscillations is increased (r, > 1) or decreased 
(r,; < 1) with respect to that for steady flows depending on the relative orders of magni- 
tude of the various terms in Eq. (9). Let us consider the first cycle (i.e., k = t). Then 
Eq. (9) can be written, after some simplification, as 


Me le a 
= | 3° vases 9 (Pr) | pra 1 + (Pr/2)'” sin (2nto 0.2233). (10) 


On the other hand, if n(k — t))/x is appreciable (over 1/2, say) then the first two terms 
of Eq. (9) give r, with an error of less than 5 per cent for air. 
For air, (Pr = 0.73) Eq. (10) shows that, for k = 0, 


0.4389 <r, < 0.593, 


depending on the value of nf) . On the other hand, for large k, the recovery factor in- 
creases as (k — f,)'’’, which increase is directly due to the energy added to the boundary 
layer by the plate oscillations. This variation is shown in Fig. 2. It is seen that advan- 
tageous recovery factors (i.e., 7; < 1) will only persist for a short time, until the initially 
cold boundary layer (Fig. 1, I) is heated by the plate oscillations. 
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A UNIQUENESS THEOREM FOR THE COUPLED THERMOELASTIC PROBLEM* 
By J. H. WEINER (Columbia University) 


1, Introduction. In the computation of thermal stresses in an elastic solid it is 
customary to compute first the temperature distribution by use of the Fourier heat 
conduction equation and then to determine the resulting thermal stresses according 
to the usual thermoelastic theory. Although this procedure is sufficiently accurate for 
a large class of problems, it is approximate since the Fourier heat conduction equation 
is an energy balance which neglects the interconvertibility of mechanical and thermal 
energy. If this possibility is included in the analysis, the energy balance equation con- 
tains both thermal and mechanical terms and the thermal and thermoelastic problems 
are coupled; the temperature and stress distributions must be determined simultaneously 
rather than consecutively. When the non-uniform temperature distribution is primarily 
due to heat supplied to the body from external sources, the mechanical coupling term 
in the energy balance may be neglected in comparison with the thermal terms; when 
the temperature differences are due solely to the deformations of the body, as in a study 
of thermoelastic damping for example, then the coupled problem must be considered. 

The coupled nature of the thermal-thermoelastic problem was already known by 
Duhamel [1]; a derivation of the governing equations utilizing thermodynamic principles’ 
was given by Voigt [2] who presents a linear theory valid for sufficiently small temperature 
changes, displacements and displacement gradients. The magnitude of permissible 
temperature changes depends upon the mean absolute temperature of the solid and the 
degree of temperature dependence of its elastic and thermal properties, while the re- 
strictions on the displacements and their gradients are those required in the linear theory 
of elasticity. 

In this note, a uniqueness theorem is presented for the coupled problem formulated 
by Voigt for the case of an isotropic elastic solid. The notation used is as follows: ¢;, , 
€;; , U; are components of the stress, strain and displacement tensors referred to a cartesian 
coordinate system z,(7,7 = 1, 2,3), A, » are Lamé’s constants, p is the density, a is the 
coefficient of thermal expansion and m = (3\ + 2y) a. T is the absolute temperature, 
T, is a reference temperature chosen so that | (7’ — T,)/T. | < 1 throughout the body, 
K is the thermal conductivity, c is the specific heat for processes with invariant strain 
tensor. The comma notation is used for derivatives with respect to space variables, 
superposed dots for derivatives with respect to the time, ¢t, 6;; is the Kronecker delta 
and the summation convention is employed. 

2. Theorem. Given a regular region® of space V + S with boundary S. Then there 


exists at most one set of single-valued functions o,;(P, f) and e,;(P, t) of class C“”, 
u,(P, t) and T(P, t) of class C” for P(z, , x2, 23) in V + S, t > O which satisfy the 


' 
following equations for P in V, t > 0, 





*Received March 8, 1956. 
1Jt should be noted here that these thermodynamic principles, as applied to deformable media, have 


not yet been put on a firm logical foundation, see [3], pp. 170-171. See also [7], which came to the author’s 
attention after submission of the manuscript, for a comprehensive treatment of this subject from the 
viewpoint of irreversible thermodynamics and a physical interpretation of the integral obtained in 
Eq. (15) below. 

As defined in [4], p. 113. 
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KT 4. = pel” + mT c€ix , (1) 
Oi3,5 = pus’, (2) 
the following equations for P in V + S, ¢ > 0, 
ei = 2(ui.5 + Uj.4), (3) 
O;; = 5,;Aex + Que;; — 6,;mT, (4) 


the following equations for P on S, ¢ > 0, 
T = FP, Ob, (5) 
u; = G;"(P, 0), (6) 


Il 


and the following equations for P in V, ¢ = 0 


T = FP), @) 
us = GP), (8) 
ui = GP), (9) 


where the constants K, c, A, u, m and T> are all positive. 
Proof. Let there be two such sets of functions, o{;’ and o{;’, e{;’ and e{;’ etc., and let 
o* = o\? — oS, e& = ef; — ef}, etc. By virtue of the linearity of the problem, it is 


clear that these difference functions will also satisfy Eqs. (1)-(4) and the homogeneous 
counterparts of Eqs. (5)-(9). In the calculations which follow the stars will be omitted 
from the designations of the difference functions. Consider the integral 


[eves dv =f osmisaV = [ Wounds — Cum av, (10) 
V “V <V 


where Eq. (3) and the symmetry of the stress tensor have been utilized. By use of the 
divergence theorem and the homogeneous form of Eq. (6), 


/ (e,;3),, dV = / o,uin; dS = 0. 
Vv s 
Also, from Eq. (2), 
a lake 0 ve > 
[ oss.su; dV = il pu;u; dV = / =, (Feuiui) dV. 
Jv v y Ot 
Therefore, Eq. (10) becomes, 


[oui + < Coun) | dV = 0. (11) 


From Eq. (4) it is found that 


0 
00:3 = NCmmCix + Qye;,e:; — mMTen = Fy (dei, + me, je:;) — mTen 








104 NOTES [Vol. XV, No. 1 
so that Eq. (11) may be rewritten in the form 
‘ re) 1 2 1 2 ry... r 
| ry, (Sreue + mesje:; + Zpu;) — mTe, | dV = 0. (12) 
JV 
The following identity is readily derived by use of the divergence theorem 


il TT xx dV + / as dV = i (TT 4). dV = [ TT im dS. 
V V V Ss 


Substitution for 7 ,, from Eq. (1) in the above and use of the homogeneous form of 
Eq. (5) yields 


[ T(pcT" + mT vis) dV + K [ T.T,dV =0 
JV JV 


or 


1. _ — ff Om _y 1 
T, [ mTei, dV = Ls? dV —K [ TT, dV. (13) 
Substitution of Eq. (13) into Eq. (12) and interchange of the order of differentiation 
and integration then yields, 


d (3 ; ae —K 
— =~ Neue + ueie:i; tu pui ta_T)dV = [ 7 cnr <6. 14 
Be te oo 27, —<_—" = alas 

The integral on the left hand side of the above equation is initially zero since the 
difference functions satisfy homogeneous initial conditions. By the inequality derived 
above, however, this integral either decreases (and therefore becomes negative) or 
remains equal to zero. Since its integral is the sum of squares, however, only the latter 


alternative is possible. That is 


[ (3 here + meses; + : nu + & 1°) dV =0, t> 0. (15) 
Jy \2 2 2T, 

It follows from Eq. (15) that the difference functions are identically zero throughout 
the body and for all time and the theorem is proved. 

3. Remarks. As noted in the Introduction, omission of the last term in Eq. (1) 
reduces the problem defined by Eqs. (1)-(9) to a heat conduction problem and a thermo- 
elastic problem which are uncoupled. In this case, the right hand side of Eq. (13) is 
equal to zero and this equation may be used directly to prove the uniqueness of the 
temperature distribution. It follows then that the difference temperature distribution, 
T, in Eq. (12) is zero and that equation may be used in the usual manner to prove the 
uniqueness of the solution of the thermoelastic problem. 

The above theorem may be generalized readily to include the most general linear 
thermal and mechanical boundary conditions’. Also the analytical restrictions on the 
solutions may be considerably lightened. These generalizations have not been con- 
sidered here since the primary purpose of this note is to indicate how methods of proving 
uniqueness of the uncoupled heat conduction and elastic problems may be combined 





3More general mechanical boundary conditions which are sufficient to prove uniqueness are, in fact, 
implicit in the second integral of the equation following Eq. (10). 
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to prove uniqueness of the coupled problem. The reader is therefore referred to other 
sources, e.g., [5, 6], for more detailed treatments of the uncoupled problems. 

Acknowledgement. The author is indebted to Professor R. D. Mindlin of Columbia 
University, who suggested the subject of this investigation. 
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THE MOTION OF A THERMOELASTIC SOLID* 
By M. LESSEN (University of Pennsylvania) 


Introduction. It is the purpose of this paper to set down formally the relevant 
equations of thermo-elasticity within the approximation of the theory of elasticity of 
infinitesimal displacements and displacement derivatives. It is not claimed that the 
following equations of thermo-elasticity are original; Duhamel [1] and Neumann [2] 
derived similar equations many years ago, but due to the fact that they did not derive 
their equations from thermodynamic considerations and also that elasticians generally 
are not aware of the role of thermodynamics in their field, it was felt that a derivation 
of the general equations and application to a particular problem were in order. The 
present work is a refinement of Refs. [3] and [4] and the application was inspired by 
the work of Synge [5] in connection with the motion of a viscous, heat conducting fluid. 
The author is indebted to the reviewer for his extensive commentary which assisted 
materially in the revised version of this paper. 

Analysis. For the case of small displacements and displacement derivatives, the 
momentum and energy (First Law) equations for a continuous, homogeneous medium 


may be written as 











Momentum p “us == Tees (1) 
at” ‘. ’ 
Energy p ; = K i oa T ns (2) 
“ at mn ~™mn mn at , 


where p is density, u, is displacement vector, ¢ is time coordinate, 7,; is stress tensor, 
U is specific internal energy, K.,,,, is thermal conductivity tensor, and 7 is the temperature. 
The subscript notation is that of cartesian tensorial form and subscripts following a 


*Received February 14, 1956; revised manuscript received May 11, 1956. 
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\ 


comma indicate partial differentiation with respect to the appropriate independent 


spatial variables. 
The momentum and energy equations constitute four equations connecting the 11 


dependent variables: 
3 te203;T 
therefore, seven additional relations must be obtained, in order to reduce the number 


of dependent variables to four. This reduction can be accomplished thermodynamically. 
For an elastic solid, where 7;;; €;;; U; and 7 are thermodynamic properties of a state, 


the Gibbs total differential equation may be written as 





T dS = qu — = oe (3) 
where S is the specific entropy, and e;; = 3(u;,; + u;,;) is equal to the “pure” strain. 
Therefore, from an assumed equation of state of the (normal) form 

U = U(S, «;), (4) 

T and 7;; can be obtained at once by 
aU oU . 
7) tea (5) 


Expanding (4) about the reference state S = S“’; e,; = 0, including second order 
terms, we obtain 


ea ) l 9 v4 | 
U = 1 ian -+ sU™ t A 55; -- 58 U™ + 8B; ;€;; + 2p Csi mn€ijEmn ’ (6) 
wheres = S— S;U, U™ ,A,; , Bi; , Csimn are the known properties 9U/ds, °U/ds’, 
dU /de;; , 0°U/de;;08, pd°U/d€; ;0€nn Tespectively at the reference state. Therefore 
T=U” +sU™ + Bi; , (7) 
T ij = pA,; aa psB;; + Crcifinn ° (8) 


U“ is seen in (7) to be the reference temperature T°’. 
Substituting (7) and (8) in (1) and (2) yields (to the first order) 


ou; 
p ry aaa pB, 58, ; + Cssenthn of (9) 
and 
TS _ KU + Buus 10 
p at _— mn\S + 6jUs.i) mn ° ( ) 


Thus, the reduction to four variables, s and u,; has been accomplished. 
The reduction to the variables JT and u; may be accomplished by solving for s in (7), 


$s = ain (T - , ae aes Biantim.n) (11) 


and substitution in (9) and (10) 
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ou,  _1 1 = os 
Ful - U® B,,T ; + (2 Css]: _— U® B.;Ban)tm.i ’ (12) 
oT — By, Mat = DasP ne, (13) 
where 
Uu 
Dan = oT | (14) 


For the case of isotropy, 
B,; = Bé;; ; |: = a ’ 
Cisne = N85; bmn + Ub in Sin + 5 in 5m) 


and (12) and (13) become 


~2 2 
=e = . 5 T 3+ (Ate - Ps inn ou, am > (15) 
Ot U p U p 

oT in 

at Fag = DT me si 


It should be noted that Eq. (16) resembles the customary heat conduction equation 
except for the term Bdu,,.,,/0t. Duhamel included a similar term in his considerations, 
but did so because he reasoned that the rate of dilatation would have a linear effect on 
the rate of temperature change. While it is true that the effect of B is small for many 
elasticity problems, many elasticians are unaware of the nature of the approximation 
they make when they use the Fourier equation to obtain a space-time temperature 
distribution. 

Propagation of waves. In the manner of Synge, let us now study solutions of 
Eqs. (15) and (16) of the form 


U; u* exp (a,x, + bt), 


T —T” = T* exp (a,x, + bt), 


where a, , 6 are complex. 
Substituting in Eqs. (15) and (16) and letting a,a, = A’; one obtains 





/ 2 
(b? a #) tus ~ (te - Br anal - ian aT* = 0, (17) 
~Bba,ut + (b — DA)T* = 0. (18) 


If a; = ai + iat’ where a} , a{’ are real, then if we choose the x; axis perpendicular to 
the vectors a{ and a’ , a, = 0. The determinantal equation for the system then is 














2_ (A+ B* d B’ B 
[ofa (te alt] -QEB— Flom  — Gers lag 
r B?’ B’ B 
-( >. ‘- Fa )avo, E = : A? — (te a ae ~pyoe a, |= 0 

—Bba, —Bba, (b — DA”) 
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or, letting 


M 72, A + 2u Pe . 72 
so a = Vr 
p p ( 
and 
A+ Qu = V2 
= V3, 
p 


the determinantal equation may be written 
(b? — V?A)[(b? — VZA?)b — (0? — V7A2)DA’] = 0. (20) 


V, V, and V,7 are clearly velocities of propagation of transverse, isentropic-longi- 
tudinal, and isothermal-longitudinal disturbances, respectively. The transverse mode is 
uncoupled from the longitudinal mode of motion. 

A few limiting cases are of interest. For B — 0; V, — V, and the mechanical mode 
uncouples from the temperature mode. The propagation of a temperature disturbance 
is then given by 


b — DA’ = 0. (21) 


For D — 0; a longitudinal disturbance is given by 


b? — VsA’* = 0. (22) 


For D — ~; a longitudinal disturbance is given by 


b? — VzA’* = 0. (23) 


The general equation for a longitudinal disturbance is seen, from Eq. (20), to be 
biquadratic in A. The two pairs of roots, therefore, represent advancing and receding 
waves of two families, each family having a different velocity of propagation. It is, 
therefore, observed that the herein studied continuum model displays a “second sound” 
phenomenon and that disturbances of both families exhibit frequency dependent damping 


and velocities of propagation. 
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A NOTE ON POISSON’S INTEGRAL* 
By R. J. DUFFIN (Carnegie Institute of Technology) 


Let D be a convex open set with boundary B. Let f be a function defined at the 
points of B and continuous at these points. This note concerns a mean value formula 
which extends f to be a continuous function in D + B. At first discussion is restricted 
to the case that D is a bounded set either in space or in the plane. 

A line drawn through a point p of D will intersect B in exactly two points, say q, 
and q, , as is illustrated in Fig. 1. The value of the function f at these points will be 





q. qo 


Fie. 1. 


denoted by f, and f, . A linear function g is defined on this line so that g takes on the 
value f, at g, and the value f, at g, . The value of g at a given point p is seen to be a 
continuous function of p and of the direction of the line through p. With p held fixed, 
let F(p) denote the average of g for all directions of the line. Thus F(p) is a continuous 
function for p in D. (It is of interest to note that if f is the boundary value of a linear 
function h, then F = h in D.) 

Let gq be a point on B where B is strictly convex. It is desired to show that if p — q@ , 
then g — f, uniformly with respect to the direction of the line through p. If this is not 
true, there is a sequence of points p, such that p, — g and | g — fo| > ¢, for a positive 
constant c, . By strict convexity | g, — qo || G2 — qo | —~ 0. Without loss of generality it 
may be assumed that g, — q . If there is a positive constant c, such that | q. — q, | > ¢2, 
then by continuity and the definition of g, it is seen that g — f, . If there is no such 
c, , then for some infinite subsequence both q, — qo and gq, — gq and so g — f, . In either 
case there is a contradiction, and so g must converge uniformly to f, . It follows immedi- 
ately that F(p) — f(qo) as p — q . This statement is seen to hold at an arbitrary point 
of the boundary by slightly modifying the above proof. 

Let r, be the distance from p to q, and let r, be the distance from p to g, . Thus 


ae 9 = (rife + refs)/( + 12). (1) 


*Received May 7, 1956. The preparation of this note was sponsored by the Office of Ordnance 
Research, U. 8. Army, Contract DA-36-061-OR D-490. 
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In the two-dimensional case 


1 2/1 
Fe) = 2 gdp = 2 et we. @) 


Here ¢ denotes the angle the line pq, makes with a fixed line. If the boundary is sufficiently 
smooth, (2) can be transformed into a line integral. Let @ be the angle formed by the 
normal to B at q, and the line pq, . Then r, dé = cos 6 ds where ds denotes the element of 


arc length of B. Thus 
1 f° cos 6 1 f° cos @ 
Fp) == [O° ds— 2 f OE fas, (3) 


where c is the length of B. It is clear that the first integral is harmonic in D. In the case 
that B is a circle of radius a, it is apparent that r, + r, = 2a cos 6. Thus the second 
integral is constant and so F is harmonic if B is a circle. 

In three dimensions formulas (2) and (3) become 











Fo) = 5 {fa an (2) 
“ii Fly) = me y. cos 6 fi ’ 
pb Sinas- thas wy 


Here dQ is the element of solid angle and dS is the element of surface area. The first 
integral in (3’) is seen to be harmonic. If B is a sphere of radius a, then again r, + r, = 2a 
cos 6 and it is seen that the second integral in (3’) is harmonic. 

Thus in the case of the circle and the sphere the mean value formula yields a harmonic 
function. Since the Dirichlet problem has a unique solution, it follows that the mean value 
formula is equivalent to Poisson’s integral. The formula may also be regarded as an 
extension of the Gauss mean value relation. There is no difficulty in extending the mean 
value formula to hold for infinite regions. This yields relations equivalent to Poisson’s 
formula for the half-plane and the half-space. 

For regions of the plane not too far from circular shape it is to be expected that 
the mean value formula (2) would give an approximate solution of the Dirichlet problem. 
Furthermore, the integral is suitable for approximation by a discrete sum. This is the 
method of the linear rosette proposed by M. M. Frocht in his book Photoelasticity, vol. 
II, John Wiley, New York, 1948. Frocht carries out some examples of the approximation 
method and makes comparison with the exact solution. 

Again consider a region of the plane. An analytic function is defined by 


we) == [ e f, de. (4) 


u=1/" ja (5) 


-! [ tan Of, dé. (6) 


Let W = U + iV then 


and 
V 


lI 
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It is seen that U is identical with the first integral in (3). Thus for a circle, U differs 
from F by a constant, and so V is the harmonic conjugate of F. 

For regions not too far from circular it is to be expected that U + C where C isa 
constant would approximate the solution of the Dirichlet problem. Of course U + C 
would not take on the correct boundary values. By the maximum principle it results 
that the greatest error would be on the boundary; hence the error could be determined. 

Since the error in the above procedure is harmonic, a second approximation could 
be set up, etc. This would lead to the method of the arithmetic mean devised by Neumann 
to solve the Dirichlet problem in convex regions. 


ON SOME EFFECTS OF VELOCITY DISTRIBUTION IN ELECTRON STREAMS 
Quarterly of Applied Mathematics, XII, 105-116 (1954) 
By S. V. YADAVALLI 


Although the details of the derivation are different Equations (83), (73) and (67) 
derived earlier by J. K. Knipp (see: “Klystrons and Microwave Triodes”, M.I.T. 
Radiation Laboratory Series, Vol. 7, McGraw-Hill, Chapter 5) are somewhat more general 
than the corresponding Equations (11), (14) and (43a) of my paper. This was not indicated 
in my above paper due to an oversight. The author would like to thank Professor J. K. 
Knipp for kindly drawing his attention to this matter. 
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BOOK REVIEWS 
(Continued from p. 82) 


Théorie des fonctions aléatoires. Applications a divers phénoménes de fluctuation. By 
A. Blane-Lapierre and Robert Fortet. Masson et Cie, Editeurs, Libraires de |’Aca- 
démie de Médecine, 120 Boulevard Saint-Germain, Paris, 1953. xvi + 693 pp. $18.77. 


The book before us sets itself the difficult task to serve the needs both of the mathematician who 
wants to know the relevant problems and of the scientist who needs and applies the theory. There is 
always a gap between theory and applications; here, however, the situation is perhaps particularly 
strained: The domain is difficult and comparatively little known: the authors rightly say that a research 
worker who may be able and willing to acquaint himself to a certain degree with a branch of analysis 
new to him will in general experience considerable difficulties and inhibitions if probability enters the 
field. (If this is valid in France, the homeland of probability calculus, it is certainly more than true in 
this country.) Random functions appear in the most varied branches of knowledge: in physics, viz. 
in the theories of turbulence, of Brownian motion, of noise, etc.; in meteorology; in economics in the 
theory of time series, and in other fields. 

The content of the 700-page work is very rich. Some of the chapters speak particularly to the 
scientist, others to the mathematician. The main topics are: Physical introduction to the theory of 
random functions (for the scientist); axioms, concepts, and theorems of probability calculus (for the 
mathematician); random functions; Poisson random process; Laplace random process; Markoff chains; 
harmonic analysis of random functions; various applications. A comprehensive very valuable chapter 
by J. Kampé de Fériet deals with the statistical theory of turbulence. The book, on the whole a great and 
successful effort, has been invited by Professor Darmois for his Series of Mathematical Works for 
Physicists. 

In this reviewer’s opinion, the book makes very difficult reading: Much is assumed as known, the 
presentation is very brief; it is interesting for re-evaluating a subject one knows well, but hard indeed 
for trying to learn new material. It is a book for the research worker who either wants to get a survey 
of existing problems and results, or who is willing to give considerable effort to the study of a particular 
subject, using other books and papers as a complement. 

The valuable bibliography contains several hundred titles,—although well known names are 


missing—; the index, hardly two pages, is meager. 
HitpaA GEIRINGER 
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VECTOR SPACES AND MATRICES 


By Rosert M. Turarr, University of Michigan; and Leonarp TornHeEIM, Cali- 
fornia Research Corporation. A new work which proceeds simultaneously on two 
levels—one concrete and the other axiomatic. The concrete approach is via matrices, 
and the intrinsic approach is via linear transformations. A unique feature of the 
book is the inclusion of a theory of algebra generated by a single element over the 
field, In the final chapter the authors give a brief introduction to game theory and 
inequalities. 1957. Approx. 320 pages. Prob. $7.00. 


ELEMENTS OF GASDYNAMICS 


By H. W. Lrepmann and A. Rosnko, both of California Institute of Technology. 
Supplies the physical background for work in modern gasdynamics and high speed 
aerodynamics. All the important fundamentals of gasdynamics are covered, in- 
cluding thermodynamics, basic equations of motion, principles of high speed aero- 
dynamics, experimental methods, viscous flow, and gas kinetics. A publication in 
the GALCIT Aeronautical Series, Theodore von Karmén and Clark B. Millikan, 
Editors. 1957. 439 pages. 167 illus. $11.00. College Edition also available. 


DYNAMIC FACTORS IN INDUSTRIAL PRODUCTIVITY 


By Seymour Metman, Columbia University. The first systematic examination of 
the causal factors of productivity. This work defines the primary causes of change, 
giving new insights into the dynamics of industrial production. It supplies explana- 
tions to differences in productivity levels of various countries and within plants in 
many industries. The author also develops methods for selecting equipment with 
an eye toward future cost trends, and for reducing managerial costs. 1956. 238 pages. 


Illus. $4.75. 


AUTOMATIC DIGITAL COMPUTERS 


By M. V. Wirkss, Cambridge University. Provides a general introduction to the 
principles underlying the design and use of digital computers. It covers the subject 
known as logical design without entering into a detailed discussion of electronic 
circuit techniques. “The reviewer . . . finds this one of the best-integrated and most 
intelligible books on automatic digital computers that has appeared to date.”—Con- 
trol Engineering. 1956. 305 pages. Illus. $7.00. 


INTRODUCTION TO OPERATIONS RESEARCH 


By C, West Cuurcuman, Russett L. Ackorr, and E. Leonarp Arnorr, all of 
Case Institute of Technology. An introductory treatment which provides a basis for 
evaluating the field and for understanding its procedures and its potential. Technical 
material has been simplified, but not at the risk of distortion. In a clear and straight- 
forward manner, the book presents a general coverage of such topics as inventory, 
linear programming, waiting line replacement, etc. 1957. 645 pages. Illus. $12.00. 
College Edition also available. 
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